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I. 


INTRODUCTION 


The  problem  of  determining  the  nature  or  properties  of  some  object  though 
a  non-invcisive  procedure  has  many  real-world  applications.  In  the  medical  field, 
a  doctor  might  want  to  determine  whether  or  not  a  patient  exhibiting  symptoms 
associated  with  a  brain  tumor  actually  hcis  such  a  tumor  without  opening  the  patient’s 
skull  to  look  inside.  In  industry,  an  inspector  might  need  to  verify  the  integrity  of 
a  pump  encased  in  a  pipe,  without  opening  the  pipe  to  do  so.  Many  other  similar 
situations  exist.  In  all  of  these  caises,  it  is  often  possible  to  determine  the  nature  of 
the  object’s  interior  by  measuring  its  density  or  some  other  physical  property,  and 
then  using  this  measured  data  to  reconstruct  the  object  in  question  in  terms  of  the 
measured  property. 

The  example  we  use  throughout  this  study  is  Computer  Assisted  Tomography 
(CAT),  the  well-known  CAT  scan.  In  addition  to  the  obvious  medical  use,  CAT  has 
important  military  applications  as  well,  such  as  jet  aircraft  engine  manufacturing. 
Here,  the  density  of  unknown  object,  say  a  jet  engine,  is  measured  by  passing  x-rays 
of  known  intensity  through  it  and  recording  the  emergent  intensities.  Given  this 
collection  of  x-ray  intensity  data,  the  density  of  the  engine  is  then  reconstructed  for 
quality  control. 

Other  medical  tomographic  applications  include  Positron  Emission  Tomogra- 
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phy  (PET)  and  Single  Photon  Emission  Tomography  (SPECT),  where  the  patient  is 
administered  a  dose  of  a  radioactive  drug,  which  collects  in  the  region  of  interest  and 
then  emits  radiation  which  can  be  detected  outside  the  patient. 

Yet  another  application  is  Ionospheric  Tomography,  in  which  the  electromag¬ 
netic  density  of  the  ionosphere  is  reconstructed.  Data  is  collected  by  transmitting 
radio  waves  from  beacons  on  the  earth’s  surface  through  the  ionosphere,  where  they 
are  detected  by  an  orbiting  satellite. 

A.  GOALS  OF  THE  RESEARCH 

Our  research  has  as  its  primary  objective  the  expansion  of  the  collection  of 
problem  types  that  can  be  approached  w'ith  a  multilevel  method.  We  use  the  image 
reconstruction  problem  as  the  vehicle  in  this  study.  Currently,  there  are  several  es- 
tal)lished  methods  for  reconstructing  the  density,  or  image  .  They  include  Fourier 
methods,  backprojection  methods,  and  algebraic  methods  -  all  of  which  have  advan¬ 
tages  and  disadvantages  in  terms  speed,  accuracy  and  scope.  We  restrict  ourselves  to 
the  algebraic  methods,  whose  major  limitation  is  that  they  are  slow,  and  attempt  to 
accelerate  them  by  incorporating  multilevel  technology.  We  accomplish  this  by  ap¬ 
plying  the  principles  of  multilevel  projection  methods  (PML)  to  the  algebraic  image 
reconstruction  problem.  A  secondary  goal  of  our  research  is  the  improved  performance 
of  an  algebraic  method.  Finally,  we  develop  a  multilevel  based  method  for  solving 
the  problem  of  Spotlight  Computed  Tomography  (CT), in  which  a  high-resolution  re¬ 
construction  is  desired  for  only  a  small  portion  of  the  image.  We  accomplish  this  by 
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applying  the  principles  of  PML  to  a  composite  grid  image  space,  developing  a  fast 
adaptive  composite  grid  (FAC)  method. 

Our  work  is  new  in  that  the  natural  pixel  discretization,  to  a  large  extent,  has 
never  been  analyzed  in  depth.  This  is  the  first  rigorous  application  of  (PML)  to  a 
problem  outside  of  partial  differential  equations,  as  well  as  the  first  application  of  FAC 
to  the  Spotlight  CT  problem,  or  of  natural  pixels  to  a  composite  grid  discretization. 
This  study  will  generally  follow  the  course  we  now  outline.  The  image  reconstruction 
problem  is  formally  posed,  and  its  relation  to  the  Radon  transform  established.  Some 
of  the  properties  of  the  Radon  transform  are  summarized,  as  are  several  inversion 
techniques.  Formally,  inverting  the  Radon  transform  solves  the  image  reconstruction 
problem. 

B.  STANDARD  APPROACH 

We  look  at  one  particular  inversion  method,  the  Algebraic  Reconstruction 
Technique  (ART),  in  greater  detail.  In  the  standard  approach  to  ART,  the  space 
containing  the  image  is  discretized  into  small  elements  called  pixels, a.nd  the  image 
density  is  assumed  to  be  constant  throughout  each  pixel.  This  approach  yields  a 
large,  sparse,  underdetermined  system  of  linear  equations,  the  solution  of  which  ap¬ 
proximates  the  image.  The  system  is  normally  solved  with  the  method  of  Kaezmarz, 
which  is  examined  and  analyzed. 
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C.  NATURAL  PIXEL  DISCRETIZATION 


Next,  we  adopt  an  alternative  discretization  based  on  natural  pixels  .  This 
discretization  was  originally  proposed  by  Buonocore  [Ref.  1],  but  a  careful  analysis 
of  the  properties  of  the  resulting  system  has  not  previously  been  performed.  This 
approach  produces  a  linear  system  that  is  square,  symmetric  and  in  general  smaller 
than  that  generated  using  square  pixels.  The  system  matrix  is  analyzed,  revealing 
a  rich  collection  of  linear  algebraic  properties.  The  rank  of  the  matrix  is  shown 
to  be  determined  by  the  x-ray  geometry  used  to  generate  it,  and  its  null  space  is 
characterized.  The  square  pixel  and  natural  pixel  discretizations  are  compared. 

D.  GAUSS-SEIDEL  ITERATION 

We  consider  the  Gauss-Seidel  iterative  method  for  solving  the  natural  pixel 
discretized  problem,  and  convergence  properties  of  Gauss-Seidel  iteration  when  ap¬ 
plied  to  this  problem  are  established.  A  spectral  analysis  of  a  typical  Gauss-Seidel 
iteration  matrix  for  this  problem  is  examined  and  serves  to  illuminate  the  numerical 
performance.  The  behavior  of  Gauss-Seidel  applied  to  several  test  systems  is  ana¬ 
lyzed,  and  numerical  results  are  presented,  along  with  several  reconstructed  images. 
This  behavior,  which  can  be  characterized  by  rapid  initial  convergence  followed  by 
stalling,  makes  the  problem  a  candidate  for  a  multilevel  approach. 
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E.  MULTILEVEL  METHODS 


A  review  of  the  traditional  multilevel  methodology  is  presented.  As  the  im¬ 
age  reconstruction  problem  is  not  traditional,  in  that  it  shares  few  characteristics 
with  problems  arising  from  elliptic  PDEs,  we  consider  more  general  multilevel  pro¬ 
jection  methods  (PML).  In  PML  [Ref.  2],  the  problem  is  discretized  by  orthogonal 
projections,  and  the  projections  themselves  implicitly  define  the  other  multilevel  com¬ 
ponents  that  make  up  the  method.  We  show  that  the  natural  pixel  discretization  is  a 
discretization  by  orthogonal  projections,  and  formally  cast  the  image  reconstruction 
problem  in  a  PML  setting.  Convergence  of  the  method  is  established,  and  its  behavior 
applied  to  several  test  systems  is  analyzed. 

F.  SPOTLIGHT  COMPUTED  TOMOGRAPHY 

Finally,  we  consider  the  problem  of  Spotlight  CT, where  a  portion  of  the  image 
is  desired  at  a  high  resolution.  Reconstructing  the  entire  image  at  high  resolution  is 
expensive,  so  a  composite  natural  pixel  discretization  at  different  levels  of  resolution 
is  developed.  The  resulting  system  matrix  is  again  analyzed  for  its  linear  algebraic 
properties.  As  in  the  case  of  uniform  discretization,  the  rank  of  the  matrix  can  be 
determined  by  the  geometry  used  to  produce  it,  and  the  null  spaee  of  the  matrix  is 
characterized.  The  composite  linear  system  of  equations  can  be  solved  using  a  block 
Gauss-Seidel  method.  This  approach  is  formally  shown  to  be  equivalent  to  the  Fast 
Adaptive  Composite  (FAC)  multilevel  method,  for  which  rigorous  theory  has  been 
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previously  developed  [Ref.  3].  Numerical  results  are  presented,  and  composite  grid 
images  are  reconstructed. 
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II. 


THE  RADON  TRANSFORM 


Consider  the  problem  of  determining  the  internal  structure  of  an  object  with¬ 
out  having  to  cut  or  otherwise  damage  the  object.  We  refer  to  such  a  problem  as 
a  reconstruction  problem,  and  it  will  be  the  bcisis  for  the  work  that  follows.  An  im¬ 
portant  category  of  reconstruction  problems  is  medical  imagery,  where  the  object  of 
interest  is  the  human  body,  or  some  particular  organ  inside  the  body.  The  profile  is 
then  used  to  reconstruct  the  object.  In  the  medical  field,  probes  includes  such  things 
as  x-rays,  sound  waves,  and  nuclear  magnetic  resonance  signals.  We  will  focus  on  the 
x-ray,  and  the  resulting  computer  assisted  tomography  problem.  For  foundational 
reading  see  [Ref.  4,  5,  6,  7,  8,  9,  10]. 

A.  COMPUTER  ASSISTED  TOMOGRAPHY  (CAT) 

If  a  mono-energetic  x-ray  is  passed  along  a  straight  line  through  some  homo¬ 
geneous  object,  then  the  intensity  of  the  x-ray  is  observed  to  decrease  according  to 
the  equation 

/  =  loe-^^ 

where  Iq  is  the  initial  intensity,  I  is  the  emergent  intensity,  and  p  is  the  linear  atten¬ 
uation  coefficient,  which  depends  on  the  material  making  up  the  object  [Ref.  11).  If 
the  x-ray  passes  through  two  different  materials,  traveling  a  distance  Xi  through  the 
first  and  a  distance  X2  through  the  second,  then  the  emerging  x-ray  will  be  attenuated 
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by 


I  = 

where  /xi  and  ^2  are  the  attenuation  coefficients  of  the  two  materials.  For  several 
materials  the  relation  is 

/  = 

For  a  non-homogeneous  object,  we  may  formally  let  the  number  of  materials  go  to 
infinity  while  the  distances  traveled  through  each  material  become  infinitesimally 
small.  Then  n  =  /i(ar),  the  linear  attenuation  function,  and  the  summation  becomes 
an  integral  over  the  x-ray  path  L,  yielding 

/  =  I^e- 

Now,  consider  paissing  many  x-rays,  along  many  paths,  through  an  object,  with 
the  x-ray  paths  directed  so  that  they  are  all  coplanar.  Then  we  may  write  the  linear 
attenuation  function  as  a  function  of  two  variables  /t(x,j/).  The  attenuation  of  one 
x-ray  is  then  given  by 

/  =  (2.1) 

where  the  line  integral  is  along  the  x-ray  path  L.  As  the  linear  attenuation  func¬ 
tion  characterizes  the  object  of  interest,  we  will  often  refer  to  it  as  the  image  to  be 
reconstructed.  Equation  (2.1)  can  be  rewritten  as 

P  =  -logiY)  =  lu{x,y)ds. 
io 

Figure  1  illustrates  the  path  of  one  such  x-ray  through  an  object. 
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Deicctor 


Source 


Figure  1.  An  x-ray  through  an  object. 


If  the  x-ray  source  and  detector  are  moved  along  parallel  straight  lines  past 
the  object  as  indicated  in  Figure  2,  then  it  is  possible  to  collect  a  set  of  attenuated 
intensities,  or  a  profile,  for  the  object  at  some  fixed  angle  4>,  as 

=  /  p{x,y)ds.  (2.2) 

L/ 


Figure  2.  A  profile  for  a  fixed  angle  (f). 
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The  bcisic  problem  of  computer  assisted  tomography  is  to  reconstruct  the  image 
^(x,  y)  from  a  collection  of  profiles  measured  at  various  angles.  Solving  this  problem 
is  equivalent  to  finding  the  inverse  of  the  Radon  transform. 

B.  THE  RADON  TRANSFORM 

Let  u(x,y)  be  an  arbitrary  function  defined  on  some  region  D  €  72.^.  If  L  is 
any  line  in  'R?,  then  the  mapping  defined  by  the  line  integral  of  u  along  all  possible 
lines  L,  also  a  function  of  two  variables,  is  the  Radon  Transform  of  u,  provided  the 
integral  exists.  Formally, 

R{L)  =  j  u{x,y)ds,  (2.3) 

where  the  domain  D  may  b<‘  all  of  72.^,  or  some  portion  thereof.  The  mapping  (2.3) 
was  first  studied  in  1917  i>y  Johann  Radon,  who  also  discovered  an  inversion  formula 
by  which  u  can  be  obtain«*d  fr«»m  R{L)  [Ref.  12]. 

Consider  a  parameterization  of  the  line  L  according  to 

p  =  X  cos  +  y  sin  <^,  (2.4) 

where  p  is  a  real  number  and  <t>  is  an  angle  measured  from  the  positive  x-axis.  Then 
(2.4)  determines  the  equation  of  a  line  through  the  xy-plane,  normal  to  the  unit  vector 
^  =  (cos  <^,  sin  ^)^,  and  a  distance  p  from  the  origin,  measured  along  Defining 
X  =  (x,y)^,  then  (2.4)  can  also  be  written  as 

p  =  x-f. 
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Using  either  [p,<f>)  or  (/», 0  ^  the  variables  of  the  Radon  transform,  (2.3)  can  be 
written  cis 

\{p,4>)  =  J^u{x,y)ds  or  =  j^u{x)ds. 

Figure  3  shows  the  geometry  of  the  Radon  transform  of  a  function  u(x,y)  in  terms  of 
p  and  4>,  where  (f)  is  the  angle  defining  a  line  normal  to  L,  the  path  of  the  x-ray,  and 
is  measured  counterclockwise  from  the  positive  x-axis.  The  parameter  p  is  the  signed 
distance  from  the  origin  to  the  line  L. 


Figure  3.  The  geometry  of  the  Radon  transform. 


It  is  sometimes  useful  to  express  the  Radon  transform,  in  terms  of  the  Dirac 
delta  function  5,  as 


[Ru]{p,<j))  =  /  u{x,y)S{p  —  X  cos  4>  —  y  sin  <f>)dxdy, 
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u(x,y)  =  ^ 


or 

lRu](p,()  =  [  u{x)5{p  -  X  •  ^)dx. 

For  some  objects,  the  Radon  transform  can  be  computed  analytically.  Consider 
a  constant  density  disk  of  radius  R  centered  at  the  origin.  Explicitly 

I,  + 

(2.5) 

0,  otherwise. 

Since  this  object  is  symmetric  with  respect  to  the  angle  (f>,  only  one  profile  is  required 
to  determine  the  Radon  transform.  Let  <f>  =  0,  so  that  the  line  integrals  are  computed 
along  lines  parallel  to  the  y-axis,  at  a  distance  p  from  the  origin.  For  values  of  \p\  >  R, 
the  lines  do  not  intersect  the  disk  and  the  profile  is  zero.  For  values  of  |^|  <  /?,  the 
transform  is 


[Ru)(/),0)  =  J 

The  symmetry  of  u(x,y)  yields 


dy  = 


1^1  <  R 


0,  otherwise. 

A  graphical  representation  of  the  Radon  transform  of  this  disk  is  given  in  Figure  4 
[Ref.  13]. 


The  Radon  transform  can  be  extended  to  higher  dimensions.  Integrating  u(x), 
for  X  €  7^”,  over  all  subspaces  of  dimension  n-1  is  also  a  Radon  transform.  For 
example,  if  n  =  3,  then  the  Radon  transform  is  the  set  of  all  integrals  of  u  over  all 
planes  in  We  will  restrict  ourselves  exclusively  to  functions  of  two  variables. 
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Radon  Translorm  of  a  Disk 


Figure  4.  The  Radon  transform  of  the  function  of  a  constant  density  disk. 

C.  PROPERTIES  OF  THE  RADON  TRANSFORM 

The  Radon  transform  operator  hats  many  properties  that  are  needed  when 
developing  inversion  technicpics.  VVe  will  outline  several  of  the  more  important  prop¬ 
erties.  A  detailed  examination  can  be  found  in  [Ref.  12]. 

1.  Linearity 

Given  two  functions  /  and  g  ,  and  two  scalars  Q  and  /?,  consider 

R{af -{■  Pg)  =  I  [of{x,y)  +  pg{x,y)]5{p-  xcos<l)-  ys\n<f))dxdy 
=  Q  I  f{x,y)6{p  — X  COS  <f>— y  Sin  4>)dxdy 

Jr^ 

+13  I  g{x,y)6{p  —  X  cos  (l>  —  y  sin  <f))dxdy 
Jr^ 

=  a\Rr\^0[Rg]. 

Thus,  the  Radon  transform  is  a  linear  operator.  This  property  is  important  in  that 
fairly  complex  objects  can  be  modeled  with  relative  ease  by  contracting  the  transforms 
of  a  collection  of  simple  objects  and  then  combining  the  results.  For  example,  to 
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analytically  compute  the  Radon  transform  of  an  annulus  centered  at  the  origin,  one 
would  compute  the  Radon  transforms  of  two  disks  of  different  radii,  and  subtract  the 
smaller  from  the  larger.  So  if 


u{x,y)  =  ^ 


1,  r|<x2  +  t/2<r? 

0,  otherwise, 
then  the  Radon  transform  of  u(x,y), hy  applying  linearity,  would  be 

^2  <  IpI  <  '  I 


[Ru\{p,4>)  - 


2\/t-2  - 

2  .  |/>|  <  r-i 

0,  tilherwise. 


Figure  5  shows  the  graphical  representation  of  the  Radtui  transform  of  the  annulus. 


ft  m  anwKfci 


Figure  5.  The  Radon  transform  of  Iht  function  of  a  constant  density  annulus. 
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2.  The  Shifting  Property 


Given  a  function  u{x),  we  consider  the  effect  of  shifting  the  argument 
of  M  by  a  vector  a.  The  Radon  transform  is 

{Ru{x  —  a)]{/>j  4>)  =  i  —  S)S{p  —  X  ■  ^)dx. 

Letting  y  =  x  —  a  results  in 

[Ru{y)]{p,(f>)  =  f  u{y)5{p  ~  {y  +  a)  -  i)dy 


Thus,  shifting  the  argument  of  u  by  a  vector  a  has  the  effect  of  shifting  the  resulting 
Radon  transform  a  distance  a  -  f  along  the  /?~axis.  The  shifting  property  allows  for 
simplified  computation  of  transforms  of  objects  not  centered  at  the  origin. 

Consider  a  disk  of  radius  r  centered  at  the  point  {a,b).  We  desire  the 
transform  of  u(T-a,y-b)  ,  where  «  is  given  by  (2.5).  Applying  the  shifting  property 
results  in 


[Ru]{p,(f>)  = 


2yJ ~  {p  ~  a  cos  (/>  —  6  sin  |p  —  a  cos  —  6  sin  ^|  <  r 

0,  otherwise. 

A  plot  of  the  Radon  transform  of  the  shifted  disk  is  shown  in  Figure  6. 


3.  The  Scaling  Property 


Consider 


[Rw](Q;p,a:|)  =  /  u{x)S{ap  —  X  '  a^)dz 
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Radon  Translonn  of  an  Shifiod  Disk 


Figure  6.  The  Radon  transform  of  the  function  of  a  constant  density  disk  shifted  a  way 
from  the  origin. 

|o| 

Tills  is  tlie  scaling  property. 

As  a  special  case  where  q  =  —1.  we  have 


[/2u](-p.-0  =  (/?u](p,0, 

so  that  the  Radon  transform  is  an  even  function  of  (p,  0-  This  eveness  is  significant 
in  a  practical  sense,  in  that  when  on  object  is  x-rayed,  only  the  angular  range  from  0 
to  TT  need  be  considered. 

The  Radon  transform  can  be  viewed  as  a  projection  operator.  For  a 
fixed  value  of  the  set  of  all  line  integrals  as  p  varies  is  a  projection  of  u  into  'R.. 
This  projection  is  identical  to  that  defined  by  the  reconstruction  problem,  as  the 
right-hand-sides  of  (2.2)  and  (2.3)  are  the  same,  hence  the  relationship  between  the 
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Radon  transform  and  the  reconstruction  problem.  If  we  define 


=  [R'^]{pA)  =  /  u{x,y)ds, 

J  Zj 

and  view  u[p,(j))  as  measured  data  obtained  from  x-raying  the  object,  then  u  can  be 
reconstructed  by  inverting  the  Radon  transform.  We  do  not  use  Radon’s  inversion 
formula  directly,  as  it  has  been  shown  to  be  numerically  sensitive  to  inaccuracies  in 
the  data  [Ref.  5],  which  in  practice  are  present  since  the  data  is  measured.  Practical 
inversion  techniques  have  been  developed,  a  few  of  which  will  now  be  overviewed. 

D.  STANDARD  INVERSION  TECHNIQUES 

There  are  several  categories  of  inversion  terhni<|ues  for  the  Kadon  transform, 
and  many  variations  within  each  category.  I'he  main  rategorie.s  are  Fourier  methods, 
backprojection  methods,  and  iterativ«*  mellunls.  This  work  will  concentrate  on  the 
latter  category,  but  we  briefly  discuss  the  i»tlier  iiietlnnl>  here. 

1.  Fourier  Methods 

The  Fourier  methods  are  ba.M*«l  on  the  relationship  between  the  Fourier 
transform  and  the  Radon  transform.  The  relationship  i.s  formalized  in  the  Central 
Slice  Theorem.  Consider  the  2-dimensional  Fourier  transform,  of  a  function 

of  two  variables, We  have 

FAD  =  r  r 

ZTT  •/— oo  J — oo 
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Next,  consider  the  one-dimensional  Fourier  transform  of  a  profile,  that  is,  the  trans¬ 
form  of  [/?u](p,  I)  taken  with  respect  to  p,  for  fixed.  This  can  be  expressed  as 


=  J1  iL 
=  ill 

=  — f  u(x)e~'^^'^dx 

2y/Tirh^  ^  ’ 

—  — 7=  /  u(x)e~‘"‘'d£ 

2V^Jtz^  ^  ' 


where  tD  =  ranges  over  all  of  'Rl.  The  Central  Slice  Theorem,  in  two  dimensions, 
says  that  the  Fourier  transform  of  a  projection  tt(p,  is  equal  to  a  constant  multiple 
of  the  two-dimensional  Fourier  transform  u{uj).  Explicitly  [Ref.  12] 


Theorem  2.1:  Let  the  image  u(x,y)  have  a  two-  dimensional  Fourier  trajis- 
form,  u{u!x,u>y),  and  a  Radon  transform,  u{p,<f>)  =  [Ru]{p,<f)).  If  Fi{Ru}  = 
u{uj,4>)  is  the  one- dimensional  Fourier  transform  (with  respect  to  p)  of  the 
profile  \Ru\{p,d>)i  then 

\/^u{Ux,Uy)  =  U{ijJ,<f>), 

where  and  d>  =  tan"*  (u;i,a;j,).  That  is,  the  Fourier  transform  of 

the  projection  of  u(x,y)  onto  the  line  in  the  direction  of  the  vector  (cos  4>,  sin  <j))^ 
is  exactly  a  slice  through  the  two-dimensional  Fourier  transform  of  u(x,y)  along 
that  direction. 


The  essence  of  Theorem  2.1  is 


\/^F2{u{x,y)}  =  Fi{Ru}  =  Fi{u{p,4>)},  (2.6) 
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where  F2  and  F]  indicate  the  two  and  one-dimensional  Fourier  transforms,  respec¬ 
tively.  Figure  7  gives  a  schematic  diagram  of  the  Central  Slice  Theorem  in  two 
dimensions.  The  interpolation  step  arises  from  numerical  implementation  using  the 
Fast  Fourier  Transfrom  (FFT). 


«(«j) 


Radon  Transform 


R«(p.«) 


2D  Fourier 
Transform 


1 D  Fourier 
Transform 


Fjiu)  (Wjj.Wy) 


taierpolaiioo 


FjtRo)  <(D) 


Figure  7.  Relationship  between  Fourier  and  Radon  transforms. 


To  invert  the  Radon  transform,  consider 

u(x,y)  =  Ff' |■^F,{u(/3.<p)}| , 

which  says  that  given  the  Radon  transform  of  sonte  image,  we  first  take  the  one- 
dimesional  Fourier  transform  of  each  profile.  This  gives  data  in  polar  coordinates. 
Fast  implementation  occurs  through  the  use  of  the  FFT,  which  requires  data  in  Carte¬ 
sian  coordinates.  Therefore,  these  results  are  interpolated  to  Cartesian  coordinates 
and  then  an  inverse  2-dimensional  Fourier  transform  is  taken  to  recover  the  image.  A 
family  of  such  Fourier  methods  exists,  based  on  how  the  interpolation  is  carried  out 
[Ref.  5,  14,  13,  15]. 
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2.  Backprojection  Methods 

Let  V’(p,  0  be  an  arbitrary  function,  where  p  =  =  x  cos  0+j/  sin  0, 

as  before.  The  backprojection  operator  B  is  defined  as 

=  2/’ 

Jo 

The  action  of  the  backprojection  operator  can  be  interpreted  as  follows.  Fix  a  point 
(x,t/).  Then  for  every  angle  0,  u(x,  y)  is  a  value  contributing  to  the  line  integral 
along  the  line  p  =  icos<^+  j/sin«/>.  That  is,  u{x,y)  is  part  of  {Ru\{p,<j))  for  every 
4>  €  [0, 27r).  Backprojection  assembles  at  (x,y)  the  sum  (integral)  of  all  values  to 
which  it  could  have  contributed,  i.e. 


^(x  cos<^  +  y  sin 


=  r  [Ru](p,4')d,l,, 

7— tt 

which  by  the  eveness  of  the  Hadoii  transform  becomes 

/^{/fu}(x,y)  =  2  /  [Ru]{p,4>)d(f). 

Jo 

Geometrically,  backprojection  is  a  form  of  image  reconstruction.  For  each  profile,  the 
values  corresponding  to  a  point  (x,y)  are  spread  over  the  image  region.  The  linear 
superposition  of  values  that  results  is  an  approximation  of  the  image.  Figure  8  shows 
the  backprojection  operation  for  two  profiles  taken  of  a  rectangular  object. 


Backprojection  by  itself  is  not  a  satisfactory  reconstruction  technique, 
as  is  evident  in  Figure  8  by  the  areas  surrounding  the  rectangle  that  have  been  shaded. 
Normally,  backprojection  is  used  as  an  intermediate  step  in  other  inversion  techniques 
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(b)  » 

Figure  8.  (a)  Two  profiles  of  a  rectangle,  (b)  Backprojection. 

that  are  very  effective.  One  of  these  techniques,  known  as  backprojection  of  filtered 
projections,  is  the  most  widely  used  reconstruction  method  [Ref.  16].  This  method 
can  be  formulated  as 

u{x,y)  = 

The  method  can  be  summarized  follows.  Given  the  Radon  transform  of  some  image, 
take  the  one-dimensional  Fourier  transform  of  each  projection  and  then  weight  the 
results  with  a  factor  |u;| .  This  weight  factor  is  the  filter.  Next,  take  the  inverse  Fourier 
transform  of  the  weighted  projections  and  then  backproject  the  results  to  recover  the 
image.  Figure  9  gives  a  schematic  diagram  of  the  backprojection  of  filtered  projectiojis 
method.  It  should  be  pointed  out  that  this  method  can  be  performed  entirely  in 
image/projection  space  without  the  use  of  Fourier  transforms,  by  using  convolutions 
instead. 
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Figure  9.  Relationships  in  backprojection  of  filtered  projections. 


This  alternate  path  to  the  image  is  also  shown  in  Figure  9.  For  ease  of 
development,  we  only  describe  the  frequency  domain  implementation. 

Another  closely  related  technique  involves  backprojecting  the  projec¬ 
tions  first,  and  then  filtering  the  backprojection.  This  method  is  known  as  the  filter 
of  backprojected  projections,  and  can  be  expressed  as 


u(x.y)  = 


Given  the  Radon  transform  of  some  image,  backproject  it  and  then  take  the  2- 
dimensional  Fourier  transform  of  the  result.  This  quantity  is  then  filtered  by  multipli¬ 
cation  by  |a;|  and  then  the  inverse  2-dimensional  Fourier  transform  is  taken,  recovering 
the  image.  Figure  10  illustrates  schematically  the  implementation  of  this  method. 


The  final  category  of  inversion  techniques  are  iterative.  They  involve  discretiz¬ 
ing  the  problem  into  a  linear  system  of  equations,  which  is  subsequently  solved  to 
recover  the  image.  We  wish  to  concentrate  exclusively  on  this  category  of  meth- 
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Radon  Transform 
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Transform 


Figure  10.  Relationships  in  filter  of  backprojected  projections. 


ods,  first  discussing  the  standard  techniques,  and  then  developing  a  more  efficient 
method.  In  the  next  chapter  we  will  develop  and  analyze  the  standard  technique  of 
discretization  by  squan  pij-t  /.s  .  as  well  as  the  iterative  method  of  Kaczmarz  to  solve 
the  resulting  linear  system. 
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III. 


ALGEBRAIC  RECONSTRUCTION 
TECHNIQUES  (ART) 

A.  DISCRETIZATION  BY  SQUARE  PIXELS 

The  algebraic  reconstruction  technique  involves  discretizing  the  Radon  trans¬ 
form  problem  Ru  =  f  into  a  system  of  linear  equations,  whose  solution  approximates 
u.  A  family  of  ART  methods  can  be  developed  bcised  upon  how  one  discretizes  the 
problem  and  solves  the  resulting  linear  system. 

The  standard  approach  is  to  discretize  the  problem  l*y  .square  pixels.  Let  u  be 
the  density  function  of  the  image  to  be  rironstruried.  ami  a.ssume  it  is  defined  in  a 
.square  region  of  unit  area.  This  unit  .squau*  is  subdivided  into  a  griil  of  smaller  ele¬ 
ments,  or  pixels  (from  picture  e/emrn/.'Land  «  is  a.ssume<l  to  be  constant  throughout 
each  pixel.  Let  the  image  be  divided  into  square  pixels  of  equal  area,  so  that  we 
are  approximating  the  continuous  solution  n  with  an  n  >  n  array  of  numbers.  These 
numbers  will  be  the  unknowns  in  the  lim-rtr  system  t»f  eipiations. 

We  define  the  geometry  at  wirn  li  x-rays  are  pa.ssed  through  the  image.  As¬ 
sume  there  is  an  array  of  N\  detectors  positioned  to  measure  the  intensity  of  the 
x-rays  after  they  have  passed  through  the  image,  and  an  array  of  A'^i  x-ray  sources 
positioned  parallel  to  the  detectors  so  that  the  path  of  an  x-ray  through  the  image 
is  perpendicular  to  the  detectors.  Assume  further  that  each  detector  measures  the 
intensity  of  only  one  x-ray,  and  that  the  sources  are  positioned  such  that  the  x-rays 
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cover  the  entire  image.  Let  the  source/ detector  arrays  be  rotated  about  the  image, 
stopping  at  M  angles  d>i,  </>2,  •  •  • ,  At  each  angle  the  array  produces  a  sampled 
profile,  yielding  a  set  of  M  such  profiles.  The  N\  detectors  on  each  of  the  M  angles 
determines  the  geometry  of  the  problem. 

Each  x-ray  passing  through  the  image  defines  an  equation  in  the  linear  system. 
Figure  11  depicts  the  x-ray  at  angle  <f>j  passing  through  the  image,  which  has  beeen 
discretized  into  an  n  x  n  array  of  square  pixels.  The  equation  generated  by  this  x-ray 
is  given  by 

n  n 

^  WIJXIJ  =  fij, 

/=!  J=\ 

where  fij  is  the  measured  intensity  of  the  x-ray,  the  x/j  are  the  unknown  values  for 
the  pixels,  and  the  wjj  are  weight  factors  which  are  non-zero  only  for  those  pixels 
through  which  the  x-ray  passed.  Observe  that 

n  n  . 

=  hi  ~  =  /  u{x,y)ds, 

/=!  J-\ 

i.e.,  the  sum  approximates  the  integral. 

There  are  N\  x-rays  at  each  of  M  angles  paissing  through  the  image,  for  a 
total  of  =  Ni  X  M  such  equations.  We  can  write  the  resulting  linear  system  as 
Kx  =  /,  where  K  :  ^  7^^. 

The  weight  factors  iw/j,  which  are  the  entries  of  K  ,  can  be  assigned  in  several 
ways,  depending  on  assumptions  made  about  the  physical  nature  of  x-rays.  First, 
assume  an  x-ray  has  no  width,  so  that  its  path  through  the  image  is  a  line.  Then  the 
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Deicctor 


weight  factors  can  be  assigned  as 


‘Wij 


1,  x-ray  peisses  through  pixel 
0,  otherwise 


We  call  this  approach  the  zero-one  discretization.  It  is  attractive  in  its  simplicity, 
but  it  has  several  drawbacks.  For  example,  if  an  x-ray  passes  through  the  center  of  a 
pixel,  or  just  through  the  tip  of  a  corner,  the  weight  factor  is  still  cissigned  a  value  of 
one,  which  intuitively  seems  inaccurate.  Also,  it  is  possible  for  a  ray  path  to  coincide 
with  the  border  of  two  pixels,  in  which  case  we  must  decide  to  either  assign  both 
pixels  a  value  of  zero  or  both  a  value  of  one.  Again,  either  decision  seems  inaccurate. 

Some  of  these  inaccuracies  can  be  overcome  by  letting  the  weight  factors  be 
defined  as  the  lengths  of  the  ray  paths  through  the  pixels.  This  approach,  which  we 
call  the  thin  ray  discretization,  corrects  the  problem  of  assigning  equal  weights  to  a 
pixel  regardless  of  whether  the  x-ray  passed  through  its  center  or  just  cut  its  corner. 
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This  approach  most  accurately  represents  the  line  integrals  of  [/?u].  However,  the 
possibility  of  x-rays  coinciding  with  pixel  boundaries  still  exists. 

Another  approach  to  assigning  the  weight  factors  can  be  developed  by  altering 
our  assumption  that  x-rays  have  no  width.  Let  the  x-rays  have  width,  so  that  a 
path  through  the  image  is  a  strip.  Then  the  weight  factors  can  be  assigned  as  the 
area  of  the  square  pixel  contained  within  the  strip.  We  call  this  approach  the  fat 
ray  discretization.  It  has  a  physical  justification  in  that  the  detectors  are  actually 
a  photgraphic  plate  subdivided  by  baffles  into  detection  bins.  Hacli  bin  will  detect 
x-rays  across  its  width,  so  the  rays  are  modeled  as  strips. 

No  matter  which  of  the  discretizations  is  used,  it  shoiilil  be  clear  from  Figure 
1 1  that  the  number  of  pixels  intersecte-d  by  an  x-ray  is  sinall  compared  to  the  total 
number  of  pixels.  Thus  the  resulting  matrix  will  in  ceneral  be  sparse.  .Also,  since  one 
would  desire  cis  high  a  resolution  image  a>  possible  without  subjecting  the  patient  to 
lethal  doses  of  radiation,  the  number  of  x  rii>*  used  is  uMially  le.ss  than  the  number  of 
pixels,  producing  a  rectangular  matrix  and  an  uiiderdet«*rmined  system  of  equations. 
A  typical  medical  application  could  invobe  n  =  .V|  =  .')12  and  M  =  180.  Thus 
the  size  of  the  problem  is  quite  large  a.s  well,  in  this  ca.se  involving  a  matrix  with  in 
excess  of  2.4  X  10*°  entries,  of  which  less  than  one  percent  are  nonzero.  Figure  12 
illustrates  the  sparsity  pattern  of  such  a  matrix. 

The  size  of  this  problem  precludes  a  direct  method  of  solution,  so  we  look  to 
iterative  methods  to  solve  the  linear  system.  As  the  matrix  is  in  general  rectangular. 
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Sparsity  Panern  ol  Kaczmarz  Matnx 


Figure  12.  ART  matrix. 


and  we  have  no  guarantee  of  non-zero  elements  on  the  diagonal,  the  classical  relaxation 
methods  such  as  Jacobi  and  Gauss-Seidel  are  not  appropriate  here.  The  method  of 
Kaczmarz  [Ref.  17]  can  be  applied  to  such  problems,  so  we  present  it  here,  and  then 
later  use  it  for  comparison. 

B.  THE  METHOD  OF  KACZMARZ 
1.  Definition  and  Properties 

Let  K  :  7?."*  define  the  system  of  equations  Kx  =  f.  Then 

given  an  initial  approximation  io  €  7^”  ,  the  method  of  Kaczmarz  corrects  the  ap¬ 
proximation  by  sequentially  adding  a  multiple  of  each  row  of  the  matrix  K  to  it.  The 
desired  multiple  is  that  which  causes  the  corresponding  component  of  the  residual, 
r  =  f  —  Kx,  to  vanish.  One  cycle  through  all  equations  is  called  a  s«;eep. Letting 
Wi  be  the  standard  b«isis  vector,  one  sweep  of  Kaczmarz  can  be  written  as 
For  i  =  1,  2,  •  •  •,  N 

Solve  (wi,  K{x  sK^Wi)  —  /)  =  0  for  s 
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Set  X  =  X  sh'^Wi 

Solving  for  s  yields 

{wi,  Kx -{■  sl\  K^Wi  —  f)  =  0, 

{wi,—f+sKl\^Wi)  =  0, 

—wj  f  +  swj  K  K^Wi  =  0, 

sw-  I\  h  Wi  =  W‘ 

- 

where  f}  is  the  component  of  the  residual  vector  r. 

We  chose  the  method  of  Kaczmarz  because  of  its  convergence  properties, 
which  simply  stated  are  that  if  a  solution  to  the  linear  system  exists,  then  Kaczmarz 
will  converge  to  it.  Formally,  we  cite  the  following  convergence  theorem  [Ref.  14]. 


Theorem  3.1:  Let  Hi  and  H2  be  real  Hilbert  spaces,  and  let  R  :  Hi  —¥  H2  be 
a  continuous  linear  operator.  Let  f  6  fh  be  given.  Assume  that  Ru  =  /  has 
a  solution.  If  uq  6  Range{R"),  then  the  sequence  generated  by  the  method 
of  Kaczmarz  converges  to  the  solution  of  minimum  norm  as  k  — >  00. 


If  /  €  Range{I\),  then  the  linear  system  Kx  =  /  satisfies  the  hypothesis  of  the 
above  theorem.  The  usual  choice  for  xq  is  the  zero  vector.  It  is  not  totally  clear  if  the 
minimum  norm  solution  is  the  best  in  terms  of  how  the  reconstructed  image  appears, 
but  because  we  can  find  this  solution  and  it  is  unique,  this  will  be  the  solution  we 
seek. 
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Geometrically,  Kaczmarz’s  method  acts  by  correcting  the  approxima¬ 
tion  in  a  direction  orthogonal  to  the  hyperplane  defined  by  the  equation  being  con¬ 
sidered.  Figure  13  shows  the  action  of  several  iterations  on  a  problem  consisting  of 
two  equations.  Observe  that  convergence  will  be  fast  when  the  hyperplanes  defined 
by  the  equations  are  nearly  orthogonal,  and  slow  if  they  are  nearly  parallel. 


Figure  13.  Geometric  interpretation  of  h’aczTnarz's  method. 

It  is  also  interesting  to  consider  a  physical  interpretation  of  the  action 
of  Kaczmarz’s  method  in  terms  of  the  reconstructed  image.  Starting  with  the  zero 
vector  as  an  initial  guess,  the  image  is  black.  The  action  of  Kaezmarz  is  to  add  a 
multiple  of  each  row  of  K  to  the  solution,  specifically  a  multiple  of  the  element  of 
the  residual  corresponding  to  that  row.  Each  row  of  K  can  be  attributed  to  an  x-ray 
passed  through  the  image,  so  the  effect  of  the  iteration  is  to  spread  a  multiple  of  that 
row  back  across  the  ray  path  through  the  image,  assigning  to  each  pixel  along  the 
path  that  portion  of  the  residual  proportional  to  that  pixel’s  contribution  to  the  ray 
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path  integral.  That  is,  the  correction  is  determined  by  the  amount  the  current 
approximation  fails  to  satisfy  the  equation,  normalized  by  the  area  of  that  strip 
integral.  Thus  the  action  of  Kaczmarz’s  method  is  a  form  of  backprojection,  which 
from  Chapter  II  we  know  to  be  a  primitive  technique  for  image  reconstruction  in  its 
own  right. 

One  sweep  of  Kaczmarz’s  method  involves  an  inner  product,  one  scalar 
multiplication  and  one  vector  addition  for  each  of  the  equations,  so  it  is  an  0(iV  x 
n^)  operation.  This  can  be  greatly  reduced  by  exploiting  the  sparseness  of  K  .  For 
example,  using  the  thin  ray  discretization  the  number  of  non-zero  entries  in  any  row 
of  the  matrix  will  not  exceed  2n\/2.  By  working  only  with  these  nonzero  entries,  the 
amount  of  work  involved  can  be  reduced  to  0{N  X  n). 

2.  Numerical  Performance 

The  method  of  Kaczmarz  is  applied  to  several  linear  systems  created 
at  different  x-ray  geometries,  and  with  assorted  right  sides  generated  both  analyt¬ 
ically  and  experimentally.  Effectiveness  is  measured  in  terms  of  the  2-norm  of  the 
residual  vector.  These  and  all  subsequent  numerical  computations  are  carried  out 
using  MATLAB  Version  4.1  running  on  a  SUN  Sparc  Station  10.  In  all  cases,  initial 
rapid  convergence  was  followed  by  stalling,  with  the  magnitude  of  the  residual  error 
well  short  of  machine  precision,  which  for  the  SUN  is  2.22  x  10~*®.  Figure  14  depicts 
graphs  of  the  norm  of  the  residual  and  the  convergence  factor  plotted  against  iter¬ 
ations.  Here  the  convergence  factor  at  each  iteration  k  is  computed  as  the  ratio  of 
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the  norm  of  the  residuals  after  the  (^■  +  sweep  to  the  residual  norm  after  the 
sweep.  Note  that  after  just  a  few  sweeps,  the  convergence  of  the  method  has  slowed 
significantly. 


Norm  of  Rotidool  vs.  Uetoltons 


Convergence  Hete 


Figure  14.  Convergence  of  Kaczmarz’s  method. 


The  geometry  of  the  example  problem  is  32  detectors  on  each  of  16 
angles,  with  a  square  pixel  discretization  using  a  grid  of  32  pixels  x  32  pixels.  Other 
geometries  were  examined,  and  this  is  a  typical  example.  The  image  being  recon¬ 
structed  is  a  brain  phantom,  constructed  by  overlaying  ellipses  and  rectangles  of 
various  sizes  and  grey  levels  inside  the  unit  square,  to  simulate  the  cross  section  of 
the  skull  and  brain.  The  data  vector  /  was  then  created  by  projecting  the  image 
through  multiplication  by  the  matrix  K.  Note  that  while  this  data  generation  is  ar- 
tihcial,  it  certainly  assures  us  that  /  €  Rangt{K)  and  that  infinitely  many  solutions 
exist.  Figure  15  shows  the  "exact”  and  reconstructed  images.  In  the  reconstruction, 
all  of  the  features  are  resolved  to  some  degree.  The  w'hite  skull  is  quite  clear,  as  are 
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the  general  shape,  intensity  and  location  of  all  the  features  contained  within. 


Figure  15.  Actual  and  reconstructed  brain  phantom  images. 


In  an  attempt  to  explain  the  performance  of  the  method,  a  thorough 
analysis  of  the  matrix  A'  is  required.  We  begin  that  analysis  with  its  singular  value 
decomposition. 

C.  ANALYSIS  OF  THE  MATRIX  K 

The  singular  value  decomposition  (SVD)  of  the  matrix  K  is 

K  =  (3.1) 

where  U  and  V  are  orthogonal  matrices,  and  S  is  a  diagonal  matrix  whose  diagonal 
entries  <t,-  are  the  singular  values  of  K  .  The  columns  of  U  and  V  are  known  cis  left 
and  right  singular  vectors,  respectively.  Singular  values  are  real  and  non-negative. 
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and  we  order  them  so  that 


>  <72  >  ■  •  •  >  (7r  >  0  •  •  •  0. 

The  number  of  nonzero  singular  values  r  <  N  equals  the  rank  of  K  .  Figure  16  shows 
the  singular  values  of  K  generated  at  32  detectors  on  each  of  20  angles,  with  the  image 
decomposed  into  a  32  x  32  array  of  square  pixels  using  the  fat  ray  discretization.  We 
will  refer  to  this  particular  geometry  as  the  standard  test  ^eomefry.Empirical  evidence 
indicates  that  this  matrix  K  is  typical. 


Figure  16.  Singular  valuts  of  a  typical  matrix  l\. 

This  figure  contains  a  characteristic  exhibited  by  the  singular  value  spectrum 
of  all  such  K  matrices  examined.  That  is,  the  spectrum  can  be  separated  into  three 
disctinct  regions,  or  bands,  as  shown  in  Figure  17.  The  first  of  these  is  the  left 
portion  of  the  spectrum,  or  resolvable  region  ,  where  the  singular  values  plot  nearly 
horizontal.  In  the  center  is  the  second,  a  narrow  region  that  we  will  call  the  near  null 
space  ,  shows  marked  decay  of  the  singular  values  that  ends  abruptly  as  the  singular 
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values  drop  off  toward  zero  in  the  third  region.  The  zero  singular  values  correspond 
to  the  null  space  of  K  ,  implying  that  the  matrix  is  rank  deficient.  In  Figure  16,  the 
resolvable  region  ranges  from  index  1  to  about  index  300,  the  near  null  space  from 
index  301  to  index  575,  and  those  singular  values  with  indices  larger  that  575  are  in 
the  null  space. 


Figure  17.  The  three  bands  of  the  singular  value  spectrum. 

Equation  (3.1)  can  be  rewritten  as 

KV  =  f/S, 

and  if  the  columns  are  equated  in  the  matrices  on  each  side  of  this  expression,  the 
collection  of  linear  systems 

Kvi  =  cr,Ui,  i  =  1  :  (3-2) 

arises,  where  Ui  and  Vi  denote  the  columns  of  U  and  V  ,  respectively.  If  Kacz- 
marz’s  method  is  applied  to  these  linear  systems  for  various  values  of  i  ,  it  is  possible 
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to  determine  which  singular  values  <t,  have  singular  vectors  Vi  that  are  slow  to  be 

reconstructed.  The  columns  of  V  are  linearly  independent  and  form  an  orthonormal 
2 

basis  for  7?.”  ,  the  space  where  the  images  live.  Likewise,  the  columns  of  U  form  an 
orthonormal  basis  for  7?.^,  which  is  projection  space.  Thus  any  right  side  can  be 
expressed  uniquely  as  and  any  solution  in  the  Range{I\)  can  be  expressed 

as  PjVj.  Then  by  solving  the  SVD  system 

Kv{  =  (3.3) 

whose  solution  is  uj,  for  each  i  ,  it  is  possible  to  determine  how  well  the  solution 
component  Vi  can  be  recovered.  The  solutions  should  be  close  approximations  to  the 
corresponding  singular  vectors  Vi.  The  quality  of  an  approximate  solution  can  be 
analyzed  by  decomposing  it  into  a  linear  combination  of  the  singular  vector  basis  as 

Vi  = 

j=i 

The  singular  vector  basis  is  orthonormal,  and  the  /3j’s  can  be  computed  as 

for  j  =  \  :  N.  (3.5) 

For  the  exact  solution  we  have 

1,  i=j 

0,  otherwise. 

A  plot  of  the  coefficients  for  the  decomposition  of  this  solution  would  be  a  spike  of 
magnitude  one  at  index  i  .  The  following  figures  are  plots  of  the  absolute  values  of 
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the  coefficients,  |y9j|,  of  approximate  solutions  after  1  and  25  sweeps  of  Kaczmarz’s 
method.  Shown  are  the  results  for  i  =  150,  a  Vj  the  resolvable  region,  and  for 

i  =  450,  a  Vj  in  the  near  null  space. 


D«oofnpoulon  e<  Appcojornalon  for  StngJat  Vakw  *150  D«€ompo«ilon  ol  Approxmaion  for  Snguiar  Valu*  9150 


Figure  18.  Plot  of  the  coefficients  of  the  decomposition  of  an  approximate  solution 
from  the  resolvable  region  after  1  and  25  iterations.  The  plots  are  magnified  to  bring 
out  detail. 


Oaoompo«ilor>  ol  ApproamMor)  for  Strtgutar  VoJuo  §450 


t>«oampo«r»on  of  Approamolon  for  Siogulor  Valuo  9450 


Figure  19.  Plot  of  the  coefficients  of  the  decomposition  of  an  approximate  solution 
from  the  near  null  space  after  1  and  25  iterations. 


37 


A  qualitative  interpretation  of  these  plots  follows.  We  observe  from  the  ex¬ 
periments  the  components  Vj  for  indices  in  the  resolvable  region  are  almost  totally 
recovered.  A  spike  of  magnitude  one  located  at  the  appropriate  index  is  clearly 
present,  along  with  a  small  amount  of  noise.  These  components  do  not  adversely 
affect  the  performance  of  the  iteration.  On  the  other  hand,  Vj  for  indices  associated 
with  the  near  null  space  the  components  are  only  partially  recovered.  There  is  no 
easily  recognizable  spike,  and  significant  noise  is  present.  These  components  represent 
the  unrecoverable  ,  or  slow  components  of  the  solution,  and  they  cause  the  iteration 
to  stall.  It  is  for  these  reasons  we  name  the  regions  resolvable  and  near  null  space.  Fi¬ 
nally,  it  should  be  noted  that  the  iteration  mixes  modes, that  is,  introduces  additional 
components  of  the  singular  value  spectrum  as  noise  into  the  approximation  wj  that 
are  not  part  of  the  exact  .solution  e,.  The  mode  mixing  occurs  in  the  near  null  space. 
This  behavior  has  serious  implications,  in  that  we  could  have  an  exact  solution  that  is 
defined  entirely  in  the  resolvable  region,  and  the  iteration  will  introduce  components 
(through  mode  mixing)  in  the  near  null  space  that  subsequently  will  be  difficult  to 
recover. 

From  our  analysis  thus  far,  if  the  solution  has  components  in  the  near  null  space 
of  K  ,  then  Kaczmarz’s  method  will  be  slow  in  recovering  them.  Also,  the  iteration 
mixes  modes,  introducing  components  in  the  near  null  speice  that  are  not  part  of  the 
actual  solution.  If  the  width  of  the  near  null  space  could  somehow  be  controlled,  we 
might  be  able  to  improve  performance.  We  first  considered  how  the  geometry  of  the 


38 


problem  influences  the  near  null  space.  It  is  determined  experimentally  that,  for  a 
fixed  number  of  detectors,  the  width  of  the  near  null  space  increases  with  the  number 
of  angles,  (as  does  the  dimension  of  the  null  space).  Figure  20  illustrates  the  singular 
value  spectrum  for  two  matrices  contructed  from  32  detectors,  using  both  4  and  16 
angles.  The  image  space  is  decomposed  into  a  32  x  32  array  of  square  pixels  using 
the  fat  ray  discretization.  For  the  4  angle  geometry,  the  near  null  space  is  defined 
by  indices  85  to  109,  while  for  the  16  angle  geometry  it  ranges  from  indices  300  to 
449.  This  is  an  increase  in  the  relative  size  of  the  near  null  space  from  21%  to  29% 
of  the  spectrum.  Intuitively,  one  would  expect  the  quality  of  the  reconstruction  to 
improve  with  a  larger  number  of  angles  (hence  more  data).  Thus  a  trade-off  likely 
exists  between  accuaracy  and  performance. 


Figure  20.  Singular  values  for  two  different  geometries. 


The  nature  of  the  Kaczmarz  method  itself  will  allow  us  to  determine  the  source 
of  the  near  null  space  components  generated  during  the  iteration.  Recall  that  during  a 
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sweep  of  Kaczmarz,  we  correct  the  approximation  by  adding  a  multiple  of  each  row  of 

the  matrix  K  to  the  approximation.  Since  the  rows  of  K  are  defined  in  the  same  space 

2 

as  is  the  image,  that  is,  'RJ^  ,  we  can  decompose  them  in  terms  of  the  singular  vector 
basis  and  discover  where  the  slow  components  originate.  An  investigation  along  these 
lines  reveals  that  strong  near  null  space  components  are  present  in  those  rows  of  A' 
corresponding  to  x-rays  that  nearly  miss  the  image  region.  It  should  be  pointed  out 
that  the  square  in  which  the  image  lives  is  the  image  region.  The  image  may  be  zero 
throughout  most  of  the  image  region,  but  the  value  of  the  image  is  immaterial.  We  are 
concerned  with  rays  that  nearly  miss  the  image  region.  A  ray  that  misses  the  image 
entirely  produces  a  row  of  zeros  in  the  matrix  and  is  not  used  in  the  reconstruction. 
The  rays  adjacent  to  these  are  those  are  of  interest  to  us.  It  is  our  conjecture  that  the 
rows  of  I\  corresponding  to  x-rays  that  nearly  miss  the  image,  i.e.  adjacent  to  rows 
of  all  zeros  corresponding  to  x-rays  that  entirely  miss  the  image,  are  the  major  source 
of  the  slow  components  in  the  near  null  space.  Figure  21  depicts  a  simple  example 
of  one  angle  of  a  fat  ray  discretization.  Rays  I  and  8  produce  rows  of  zeros  as  they 
entirely  miss  the  image.  Rays  2  and  7  are  near  misses  and  most  likely  correspond  to 
rows  in  the  matrix  which  contain  near  null  space  components.  The  remaining  rays 
generate  rows  that  most  likely  contain  only  components  in  the  resolvable  region. 

We  now  return  to  the  standard  tets  geometry  of  32  detectors  and  20  angles. 
Figure  22  depicts  the  spectral  decomposition  of  a  row  corresponding  to  an  x-ray 
passing  through  the  center  of  the  image  (row  8),  as  well  as  that  of  a  row  corresponding 
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Detector 


Source 


Figure  21.  One  angle  of  a  fat  ray  discretization. 


to  a  near  miss  (row  28),  taken  from  a  matrix  with  geometry  32  detectors,  4  angles, 
and  32  x  32  square  pixels.  Recall  from  Figure  20  above  that  for  this  geometry  the  near 
null  space  ranges  from  indire.s  85  to  109.  The  presence  of  near  null  space  components 
is  obvious  in  the  second  plot.  One  solution  to  this  problem  is  to  drop  the  near  miss 
rows  from  the  matrix,  but  <loing  so  results  in  a  geometry  that  does  not  completely 
cover  the  image.  Another  approach  might  be  to  keep  the  rows  in  /\ ,  but  not  include 
them  in  the  Kaezmarz  sweep.  Another  is  to  join  near-miss  strips  with  their  interior 
neighbors.  In  all  cases,  inaccuracies  might  arise  if  significant  parts  of  the  image  lie 
along  its  edges  or  in  its  corners. 


Finally,  we  consider  the  spectral  decomposition  of  the  image  we  are  trying 
to  reconstruct.  Let  u(x,y)  be  the  image  to  be  reconstructed,  and  assume  we  have 
discretized  it  into  an  n  x  n  array  of  square  pixels  x  6  7^”  so  that  the  value  in 
each  square  pixel  is  constant.  This  array  is  the  exact  solution  that  we  are  trying 
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SptcOli  DfcompoMWi  «l  Row  2ft 


Figure  22.  S-pectral  dtcomjiosiiion  of  an  inierior  and  mar  miss  row  of  K. 


to  approximate.  Based  on  experimental  data,  if  we  dm >111  pose  x  in  terms  of  the 
singular  vector  basis  of  A’, we  find  tliat  the  exart  sol  11  lion  contain.s  romponents  in 
NS({\)  .  Figure  23  shows  the  spectrum  of  a  piiantuni  brain  image  using  the  SVD  of 
a  matrix  vvith  geometry  32  detectors.  ^  aimles  and  32  •  .32  pixels.  'I'lie  null  space  of 
ill  is  matrix  is  a.ssociated  with  indice.*'  22ti  to  2'ib  h  tlear  from  the  figure  that  the 
exact  solution  has  components  in  the  mdl  vp.ue  of  the  matrix.  U'e  know  that  the  data 
vector  /  for  an  image  is  in  the  Aarir/fi  A  1  bv  tlie  way  it  is  generated,  i.e.  projecting 
it  with  A',  so  inconsistency  is  not  the  sourte  nf  tlie  problem.  These  components  in 
NS(K)  cannot  be  recovered  by  the  iteration.  they  form  a  portion  of  the  error  in 
the  approximation  that,  essentially,  ran  iii>t  be  rt'solved. 

We  may  summarize  our  analysi.s  of  the  matrix  A'  resulting  from  a  square  pixel 
discretization  of  the  image  reconstruction  problem,  noting  that  the  analysis  has  re¬ 
vealed  some  significant  drawbacks.  First,  the  discretization  results  in  a  large,  rect- 
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angular,  sparse,  singular  matrix  to  which  the  classical  iterative  methods  cannot  be 
applied,  so  we  use  the  method  of  Kaczmarz  as  our  solver,  which  tends  to  stall  after  a 
few  iterations.  Further  analysis  using  the  SVD  of  K  reveals  the  spectrum  of  the  ma¬ 
trix  can  be  separated  into  three  bands,  the  resolvable  region,  near  null  space,  and  null 
space.  We  find  that  Kacamarz’s  method  cannot  ecisily  resolve  components  in  the  near 
null  space,  and  that  it  mixes  modes,  thereby  introducing  components  in  the  near  null 
space  that  subsequently  cannot  be  recovered.  It  is  also  found  that  the  exact  solution 
itself,  when  decomposed  in  terms  of  the  singular  vector  bzisis,  hzis  components  in  the 
null  space  that  can  not  be  recovered  by  the  iteration.  In  the  next  chapter  we  develop 
a  different  discretization  of  the  problem  that  may  alleviate  some  of  these  problem 
areas. 
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IV.  NATURAL  PIXEL  DISCRETIZATION 


A.  DERIVATION 

Let  u(x,y)  be  the  density  function  of  the  image  to  be  reconstructed,  and  assume 
it  is  defined  in  a  square  region  of  unit  area.  Assume  we  have  an  array  of  detectors  to 
capture  the  intensity  of  the  x-rays  after  they  have  passed  through  the  image.  We  do 
not  require  the  detectors  to  be  evenly  spaced,  but  in  general  they  are.  Further,  assume 
the  x-ray  sources  are  arrayed  parallel  to  the  detectors,  so  that  the  x-ray  paths  are 
perpendicular  to  the  detectors.  The  sources  are  not  point  sources  and  the  detectors 
are  not  point  detectors,  but  have  non-zero  widths,  so  that  an  x-ray  passing  through 
the  image  defines  a  strip  and  its  emerging  intensity  is  entirely  detected  by  the  ray’s 
corresponding  detector.  Finally,  assume  that  the  sources  and  detectors  are  positioned 
so  that  there  is  always  total  x-ray  coverage  of  the  image  by  the  strips.  Now,  let  the 
.source-detector  arrays  be  rotated  through  a  set  of  M  angles,  and  a  profile  measured 
at  each  angle.  We  again  hope  to  reconstruct  the  image  density  u  from  this  collection 
of  M  profiles. 

Next,  assume  that  all  sources  and  detectors  have  a  constant  width  so  that 
each  x-ray  passing  through  the  image  defines  a  strip  of  constant  width  from  source  to 
detector,  as  depicted  in  Figure  24.  The  collection  of  strips  at  a  given  angle  completely 
covers  the  image,  and  can  be  thought  of  as  a  set  of  pixels  for  that  angle  which  are 
uniquely  defined  by  the  x-ray  paths.  For  the  M  angles  we  have  a  collection  of  these 
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pixels  that  are  discrete  and  overlapping.  As  they  arise  naturally  as  a  result  of  the 
geometry  used  to  x-ray  the  image,  they  are  often  refered  to  as  natural  pixels  [Ref.  1]. 


Figure  24.  X~ray  strips  complclt  hj  rort  ring  tin  nmigt  at  a  jiitd  angle. 

Discretization  of  the  problem  rei|uir«*>  the  introdiirtion  of  rliaracteristic  strip 
functions  corresponding  to  the  natural  pixels.  Define  .Vijm),  for  m  =  1  :  M,  to  be 
the  number  of  sources/detectors  for  the  rn"’  profile  whose  x-rays  pass  throught  the 
image.  Then  N  =  ^^=1  ^i(^)  *'•••‘1  number  c»f  x-rays  passing  through  the 

image,  and  hence  the  total  number  of  natural  pixels,  for  a  particular  geometry.  Let 
r/fk  :  1  <  k  <  N,  be  the  k‘^  chararterislir  strip  function,  where 

I  1  for  {x,i/)  €  5*. 

V’it  =  S 

0  elsewhere 
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and  Sk  is  the  region  of  the  strip  within  the  image  square.  Figure  25  depicts  one 
such  characteristic  strip  function  for  a  given  angle  (j>.  Note  that  the  strip  function  ipk 
is  nonzero  only  on  the  shaded  portion,  not  the  entire  strip. 


Figure  25.  A  representative  characteristic  strip  function. 


Define  the  operator  A  :  H  by 


Au 


\ 

(V’a.ti) 


Here,  (•,  •)  is  the  standard  L2  inner  product  and  H  is  some  appropriate  Hilbert  space 
in  which  the  image  u  is  defined. 
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The  detectors  measure  x-ray  intensity  after  the  ray  passes  through  the  image, 


producing  a  data  vector  /,  which  can  be  modeled  as 


Au  = 


N 

\ 

/l 

(V’2,u) 

= 

/a 

Jn  j 

The  system  Au  =  /  is  underdetermined,  so  that  if  it  is  consistent  then  there 
exist  infinitely  many  solutions.  We  must  select  one  solution  for  the  image,  so  we 
choose  the  minimum  norm  solution.  This  is  given  by  A’a  (Ref.  18],  where  a  € 
solves  the  system 

AA‘a  =  /. 

Once  we  have  found  the  vector  a,  then  the  minimum  norm  solution  is  given  by 


u(j’,y)  =  A’a. 


It  is  eaisy  to  see  that  the  adjointoperator  A’  :  —¥  H  is  given  by 


\ 

0^1 


A’a  =  •  •  •  M 


02 


N 

k=\ 


ii(x,  y). 


Hence,  A’  can  be  viewed  as  a  backprojection  operator.  It  cissigns  a  value  Oj  to  each 
strip  function  xf^j  in  the  image.  The  strips  overlap  and  the  Oj  values  acculumate 
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additively  in  the  intersections  of  the  strips,  ultimately  producing  a  representation 
of  the  image.  The  image  density  u  represented  as  a  linear  combination  of  these 
characteristic  functions  by 

N 

H  y)  (4.1) 

k=:l 

defines  a  grid  of  polygons,  on  each  of  which  u  is  constant.  A  typical  grid  of  these 
polygons  is  depicted  in  Figure  26,  which  results  from  a  geometry  of  20  angles  and  32 
detectors  per  view. 


Polygon  Gnd 


Figure  26.  A  representative  grid  of  polygons. 


Define  the  operator  B  :  'R.^  — f  'R,^  as 

B  =  AA\ 


Then  finding  the  minimum  norm  solution  given  by  A*a  is  equivalent  to  solving  the 
linear  system  Ba  =  /.  It  should  be  noted  that,  unlike  the  case  of  overdetermined 
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systems,  this  formulation  does  not  square  the  condition  number  of  the  operator  [Ref. 
19]. 

The  entries  in  B  can  be  calculated  by  substituting  (4.1)  for  u  in  Au  =  /, 
yielding 


N 

=  /i,  i  =  1,  ••• 

it=l 

Hence,  the  (j,  ky^  entry  of  B  is  given  by 

Tlie  problem  of  reconstructing  an  image  from  samples  of  its  profiles  is  now  discretized 
hy  Bo  =  /,  where  R  =  and  the  approximate  solution  to =  f  'lsu  =  A'a. 

The  following  simple  example  will  serve  to  illustrate  the  discretization  process. 
.Assume  we  have  two  profiles,  one  with  two  detectors  and  the  second  with  three 
detectors,  oriented  and  numbered  as  in  the  figure  below.  Note  the  the  image  space  is 
of  unit  area,  so  that  the  areas  of  1  and  2  are  ^  and  the  areas  of  4,5  and  6  are  5. 

The  matrix  B  generated  from  this  geometry  will  have  entries  corresponding  to 
the  areas  of  intersection  of  the  strips,  e.g.  element  (/?23)  will  be  the  area  of  intersection 
of  strip  2  with  strip  3,  which  is  |.  The  complete  matrix  is 
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View  1  View  2 

Figure  27.  Geometry  for  two  profiles. 
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We  must  only  solve  the  linear  system  resulting  from  this  process  for  a  to  reconstruct 
the  image  u. 

Recall  that  the  continuum  problem  we  are  trying  to  solve  is  Ru  =  f.  In  our 
discretization,  we  first  approximate  R  with  A  ,  yielding  Au  =  f .  We  use  a  finite 
number  of  strip  integrals  taken  at  a  finite  number  of  angles  around  the  image  as  an 
approximation  to  the  Radon  transform.  Next,  u  is  expanded  as  a  linear  combination 
of  characteristic  strip  functions,  yielding  AA*a  =  Ba  =  /.  Thus  a  linkage  exists 
between  the  continuum  problem  and  our  discretization. 
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We  examine  the  strength  of  this  linkage  by  examining  the  grid  of  polygons 
produced  by  superimposing  the  strip  functions  as  depicted  in  Figure  26.  Define  Cl 
as  the  unit  square  containing  the  image,  and  let  g{x,y)  be  a  continuous  function 
defined  on  Cl.  Let  C  be  the  total  number  of  polygons  in  the  partition  of  Cl,  generated 
by  strips  of  constant  width  over  equally  spaced  angles,  and  let  Ci  be  the  area  of  the 
polygon,  so  that  a;  =  1.  Finally,  let  ^{{x,y)  be  a  characteristic  function  for 
the  polygon,  so  that 


1,  (x,j/)€  z"*  polygon 

0,  otherwise. 


Then  we  have  the  following  result,  partially  attributed  to  Rhoden  [Ref.  20]. 


Theorem  4.1:  For  any  g(x,y)  and  any  c  >  0  there  exists  a  strip  discretization 
using  (  strips  and  a  function  f  =  cti^i{x,y)  such  that  \\f  —  ^|1/?  <  c, 

where 

Proof:  Since  g(x,y)  is  continuous,  there  exists  a  S  for  every  £  such  that 
||(2;i,yi)  -  (^2,J/2)||  <  ^  implies  |y(x,,t/,)  -  g{x2,y2)\  <  Let  c  be  given,  and 
let  So  be  the  S  that  implies  the  continuity  conditiion.  Let  the  strip  width  used 
in  the  discretization  be  which  in  turn  determines  the  number  of  polygons, 
which  we  denote  (q.  Then  for  two  angles  (which  is  the  fewest  allowable)  the 
grid  consists  of  squares  with  a  maximum  chord  length  of  ^o-  Adding  additional 
angles  cannot  increase  the  maximum  chord  length.  Consider 


\\f-g\\%  = 
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Since  g(x,y)  is  continuous,  it  attains  a  maximum  Mi  and  a  minimum  m,  on 
each  polygon.  Choose 

Mi  +  m, 


a,  = 


Then,  by  the  intermediate  value  theorem  there  exists  a  point  {Xi,y,)  in  each 
polygon  so  that  =  a,-.  Therefore  we  have 


(x,y)  -  (x.,y.)||  <  (Jo  \gix,y)  -  o,l  <  t. 


So 


Wf-gWn  = 


< 


I 


Hence  we  can  approximate  any  contirnu»ii«>  fuiirtion  arliitrarily  closely  on  a  polygon 
grid  generated  by  strip  pixels.  Since  the  rontimioiis  functions  are  dense  in  L2  [Ref.  21], 
this  result  implies  that  we  can  arbitrarily  closely  approximate  any  function  defined  on 
^2(0)  by  any  continuous  function,  which  in  turn  can  be  approximated  by  a  function 
defined  on  the  polygon  grid  to  arbitrary  precision. 

We  hope  to  find  that  the  linear  system  produced  by  discretizing  the  problem 
using  natural  pixels  will  have  distinct  advantages  over  that  produced  using  a  square 
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pixel  discretization.  Specifically,  we  seek  a  system  that  lends  itself  to  other  itera¬ 
tive  solvers  than  Kaczmarz’s  method,  one  that  requires  less  work  to  solve  while  still 
yielding  an  image  of  comparable  quality,  and  one  that  might  be  approached  with 
a  multilevel  method.  A  careful  analysis  of  the  matrix  B  will  help  us  realize  these 
advantages. 


B.  ANALYSIS  OF  THE  MATRIX  B 

The  matrix  5  is  N  x  N  with  entries  {fijk)  given  by  for  j,k  =  I  :  N. 

The  quantity  is  just  the  area  of  intersection  of  the  and  x-ray  strips. 

Figure  28  illustrates  the  geometry  that  would  generate  entry  (/3jfc)  of  the  matrix. 


Area  =  (bj|() 


Figure  28.  Intersection  of  two  x-rays  inside  the  image. 


The  nature  of  the  non-zero  entries  of  B  gives  the  matrix  a  rich  collection  of 
properties  that  we  can  exploit,  the  first  of  which  is  symmetry. 


Lemma  4.1:  B  is  symmetric  and  nonnegative. 
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Proof:  B  =  (/?,j)  =  {xpi,rpj)  =  =  (/3j«)  =  The  non-zero 

entries  of  B  represent  areas  of  polygons,  and  as  such  cannot  be  negative.  | 


If  B  were  positive  definite,  then  any  of  the  classical  iterative  methods  could  be 
applied  to  the  linear  system,  with  convergence  guaranteed  [Ref.  22].  We  can  show 
that  B  is  positive  semidefinite,  and  later  use  this  fact  to  show  that  the  Gauss-Seidel 
method  when  applied  to  the  problem  cannot  diverge  and  in  fact  must  converge. 


Theorem  4.2:  B  is  positive  semidefinite. 

Proof:  Let  x  €  7^^  be  any  non-zero  vector.  Then  we  can  write 

N  N  N 

x'^Bx  =  XiY,  +  a:2  XI  -|-  •  •  •  +  X 


1  =  1 


i=l 


1=1 


N  N  N  N 

J=1  1=1  J=1  1=1 


But  0ji  =  {xl)j,^i),  so 


N  N 

Bx  = 

i=i 1=1 
N  N 

j=\  1=1 


/  N  N  \ 

\j=l  i=l  / 

Therefore,  B  is  positive  semi-definite. 


N 


X 


1=1 


>  0. 


B  has  a  special  block  structure,  and  can  be  expressed  as 
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\ 

•••  Bim 
B22  '  •  •  B2M 

J 

^  Bm\  Bm2  •••  Bmm  ^ 

where  each  block  results  from  the  intersections  of  all  the  rays  within  two  given  angles. 
If  we  define  the  angles  used  for  the  M  profiles  to  be  4>\,02^  4>Mi  then  block  Bij 
is  formed  by  considering  the  intersections  of  the  rays  at  angle  o,  with  those  at  angle 
<i>j.  The  size  of  the  block  is  the  number  of  detectors  at  angU*  o,  by  the  number  of 
detectors  at  angle  4>j,  or  A^j(f)  x  Ni{j).  F'igure  25)  illustrates  the  block  structure  of  a 
typical  matrix,  resulting  from  a  geometry  of  S  angU*s  and  «h*tectors  per  view. 


B  = 


Bn 

B21 


Block  Strucfurv  if  mrc  • 


ru  *6662 

Figure  29.  The  block  structure  of  B.  Sou-zero  elements  are  highlighted. 


When  angle  equals  angle  <^j,  the  only  time  that  the  two  rays  will  intersect 
is  when  they  coincide.  Therefore,  the  blocks  Bu  are  diagonal  matrices  whose  entries 
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are  the  areas  of  the  natural  pixels  at  angle  4>i-  The  natural  pixels  from  four  detectors 
at  a  45  degree  angle  are  shown  in  Figure  30. 


Figure  30.  Natural  pixels  from  four  detectors  at  angle  <f>i  =  ^5  degrees. 

The  resulting  diagonal  matrix  is 

/  \ 

I  0  0  0 

0^00 

Bii  = 

0  0  I  0 

^0001^ 


Lemma  4.2:  The  elements  of  the  blocks  of  B  exhibit  the  following  summability 
properties: 

a)  The  elements  of  any  diagonal  block  Bn  sum  to  the  area  of  the  image. 

b)  Let  rit  =  \(dk\  (^k2  •  •  •  PkN]  the  k*^  row  of  B,  and  suppose  that  the 
strip  pixel  occurs  within  the  set  of  strips  at  angle  The  elements  of  rk  in 
any  off-diagonal  block  Bij  sum  to  bkk,  the  value  of  the  diagonal  element  of 
rk- 
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c)  Let  Ck  =  [/?iA:  /52A;  ' ' '  ^Nk]  be  the  column  of  B,  and  suppose  that  the 
strip  pixel  occurs  within  the  set  of  strips  at  angle  (f)j.  The  elements  of  Ck 
in  any  off-diagonal  block  Bij  sum  to  hkk,  the  value  of  the  diagonal  element 
of  Ck. 

Proof:  At  angle  <f)i,  diagonal  block  Bn  has  entries  (3kk  =  (i’k^i’k)  = 

fg‘tpl(x,y)dxdy  =  fg^dxdy.  Thus  ^kk  is  the  area  of  the  strip  correspond¬ 
ing  to  V’A:-  By  definition,  the  rays  for  any  angle  completely  cover  the  image. 
Combining  these  facts  proves  part  a  .  To  prove  part  b  ,  consider  how  the 
elements  of  row  k  for  the  off-diagonal  block  Bij  corresponding  to  angle  (f>i  are 
constructed.  The  entries  of  row  k  for  this  block  are  the  areas  of  intersection 
of  strip  k  at  angle  with  all  the  strips  for  angle  (j)j.  We  know  that  for  angle 
(f)j  the  rays  must  completely  cover  the  image,  so  they  must  completely  cover 
strip  k  as  well.  This  geometry  is  illustrated  in  Figure  31.  Therefore,  the  sum 
of  these  intersections,  which  are  the  elements  of  row  k  ,  must  sum  to  the  area 
of  ray  k  ,  which  is  the  diagonal  element  of  row  k  .  Since  B  is  symmetric,  part 
c  follows.  I 


Figure  31.  Geometric  interpretation  of  summability  property. 


We  know  from  Theorem  4.2  that  B  is  positive  semidefinite,  so  it  can  have  zero 
eigenvalues.  The  summability  properties  discussed  above  give  us  an  insight  as  to  the 
number  of  zero  eigenvalues.  Consider  any  row  of  blocks  Bn.,  Bi-i  « •  • ,  BiM,  of  5, and 
sum  the  indivdual  rows  into  a  new  composite  row.  Since  the  elements  of  the  k^^  row 
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or  column  of  any  off-diagonal  block  Bij  sum  to  the  corresponding  diagonal  element 
bkk  of  that  row  or  column,  this  composite  row  will  consist  of  the  N  diagonal  elements 


of  5  .  As  this  property  holds  for  all  M  rows  of  blocks,  B  has  at  least  a  rank  deficiency 
of  M-1  .  We  will  ultimately  show  that  the  rank  deficiency  is  exactly  M-1, hut  first 
further  analysis  of  the  structure  and  properties  of  B  is  required. 

Let  B  €  be  the  block  matrix  with  square  diagonal  blocks  of  dimensions 

ni,  7^2)  and  so  that  W]  ^2  +  ”3  =  N,  given  by 


B  = 


\ 

Bn  Bi2  Biz 


B21  B22  B23 


Bzi  B32  Bzz 


and  let  v  £  'R-^  be  given  by 


V  = 


ai 

«2 


03  / 


where 


1 


=  ai 


\ 

1 


,  02  =  02 


\ 

1 


,  and  0*3  =  03 


and  dimensions  of  the  vectors  a,-  correspond  to  the  dimensions  of  the  blocks  of  B,i.e. 
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the  length  of  a,-  is  n,.  Then 


/  ^  ^  \ 
(5iidi  +  5i2<i*2  +  B\z(^z) 


Bv  = 


(521  <5*1  +  522<i*2  +  523<53) 

^  (531<5i  +  532(5*2  +  533(53)  ^ 

We  wish  to  define  vectors  with  the  properties  of  v  in  the  above  example  as 
being  constant  over  each  of  the  M  angles  used  in  the  generation  of  the  block  matrix 
5  .  Formally: 

Definition  4.1:  A  vector  V  £  is  constant  by  angle  with  respect  to  the 
block  matrix  5  €  generated  over  M  angles  if  it  consists  of  M  subvectors 

of  constant  value  <a:i,  02,  •  •  • ,  ocm  corresponding  to  the  block  structure  of  5  . 

Consider  a  vector  v  G  that  is  constant  by  angle.  If  v  represents  the  set 
of  coefficients  defining  a  reconstructed  image  in  terms  of  its  natural  pixels,  then  for 
each  of  the  M  angles  <f>i  all  of  the  strips  covering  the  image  have  the  same  value  a,. 
Thus  for  each  angle  the  contribution  to  the  image  is  constant,  and  hence  the  image 
over  all  the  angles  is  constant,  with  a  value  equal  to  the  sum  of  the  q,’s. 

Writing  5 


5  = 


\ 

/?11 

•  •  •  filN 

/?21 

/?22 

•  •  •  fi2N 

= 

hi 

y  Pni 

Pn2 

•  •  •  ^NN  j 

we  have  the  following  lemma: 
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Lemma  4.3:  If  v  E  is  constant  by  angle  with  respect  to  the  block  matrix 

B  generated  over  M  angles,  then  bi  v  =  +  0:2  +  ■■■  +  cum),  where  the 

Q,-  are  the  constants  in  the  definition  of  constant  by  angle.  If  z  ^  is  the 

-T  _ 

vector  of  all  ones,  then  bi  z  =  [M)l3ii. 

Proof: 

-T  _ 

bi  V  =  OiAi  +  +  •  •  •  +  ai/5t;Vi(i) 

+02/3tWi(l)+l  +  <^2PtNi{l)+2  H - 1-  Ol2PiNi{l)+Ni[2)  +  '  '  ' 

+OlMPiN-Ni{M)+l  +  OlMPiN-Ni{M)+2  H - +  OlM^iN 

Ni(l)  Ni{2)  N 

=  oi  ^  fdik  +  0L2  (dik  +  •  •  •  +  oiM  E 

k=l  k=Ni{l)+l  k=N-Ni{M)+l 

which  after  applying  Lemma  4.2b  can  be  written  as 

-.T_ 

bi  V  =  aiPii  + - h  QA/Ai 

=  /3ii{oil  +  02  +  •  •  •  +  cum)- 

Since  z  is  constant  by  angle  with  all  the  o,-  =  1,  the  result 

bi^  z  =  (M)Ai 

follows  immediately.  | 

Armed  with  Definition  4.1  and  Lemma  4.3,  we  are  prepared  to  characterize 
the  null  space  of  B.  Indeed,  we  can  show  that  vectors  in  the  null  space  NS(B)  are 
constant  by  angle  and  thus  correspond  to  constant  images.  We  accomplish  this  with 
the  following  theorem,  a  portion  of  which  is  due  to  Limber  [Ref.  23]. 

Theorem  4.3:  v  €  NS{B)  if  and  only  ifv  is  constant  by  angle  and  — 

0. 

Proof:  Let  v  be  constant  by  angle,  and  let  =  0-  Then,  by  Lemma 

4.3, 

/  \ 

B5=  =0. 

V  f^NN  IZjLl  / 
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Therefore  v  €  NS{B). 

Now,  let  V  €  NS{B).  We  first  show  that  NS{B)  =  NS{AA*)  =  NS{A*). 
So,  assume  v  €  NS{A*).  Then  v  €  NS{AA*)  and  AA*v  =  0.  If  u  €  NS{A*), 
then  V  €  NS{AA*)  since  A{A*v)  =  j4(0)  =  0.  There  are  no  u  ^  NS{A*) 
such  that  AA*v  =  0,  because  if  AA*v  =  0,  then  A*v  €  Range{A*).  We 
know  Range{A*)  ±  NS{A)  and  therefore  A*v  ^  NS{A)).  Then 


A*v  =  [ipi{x,y)  ■ip2{x,y)  •••  V’jv(a;,y)] 


/  \ 
V2 

\^N  / 


N 


=  =  0. 

j=i 


Consider  two  adjacent  strips  Sk  and  Sk+i  from  the  same  profile.  Now,  select 
points  (a:i,yi)  €  D  and  (2:2, 1/2)  €  Sk+i  H  Q,  where  fi  is  the  intersection  of 
(M  —  1)  strips,  one  from  each  profile  and  none  from  the  profile  containing  Sk 
and  Sk+i’  This  selection  can  always  be  done.  To  see  why,  superimpose  all  of 
the  x-ray  path  strips  over  the  image  at  once.  They  subdivide  the  image  into  a 
collection  of  polygons.  Since  each  profile  completely  covers  the  image,  a  point 
in  the  interior  of  any  polygon  is  contained  in  one  strip  from  each  of  the  M 
profiles.  A  point  on  any  edge  (not  a  vertex)  of  any  polygon  will  be  contained 
in  Cl,  with  the  edge  separating  strips  Sk  and  Sk+i-  Moving  a  distance  e  to 
either  side  and  perpendicular  to  the  edge  will  define  the  points  {xi,yi)  and 
{x2,y2)‘  Figure  32  illustrates  this  geometrically. 


Figure  32.  Geometric  representation  of  the  proof  of  Theorem  4-3. 
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Since  V  G  NS{A'‘),  we  have 


N  N 

Y^Vi'tl>i{xi,yi)  =  Ylvii’i{x2,y2)  =  0. 

t=l  «=1 

The  ^i(x^y)  are  characteristic  functions,  so  we  can  rewrite  the  above  sums  as 

Y,vj  +  vk  =  Y^vj  +  vk+i  =  0, 

jeii  ien 

which  implies  Vk  =  Ufc+i-  Therfore,  v  is  constant  by  angle.  Further,  since 
Sjen  ‘k’j  +  Vfc  =  0  and  there  is  a  Vj  from  each  of  the  M  profiles  in  the  sum, 
the  constants  from  the  definition  of  constant  by  angle  sum  to  zero.  | 


An  immediate  consequence  of  Theorem  4.3  is  that  we  can  show  the  dimension 
of  NS{B)  is  exactly  M-1. 


Theorem  4.4:  Let  B  €  be  the  block  matrix  generated  over  M  angles 

as  discussed  above.  Then  the  dimension  of  NS{B)  =  M  — 

Proof:  Let  v  G  NS{B).  Then  v  is  constant  by  angle  and  oci  —  0-  We 
can  select  Oi,  Q2,  •••,  Qm-i  arbitrarily,  and  then  qm  is  determined.  Hence 
we  have  Af  —  1  degrees  of  freedom  in  the  selection  of  the  a,’s,  which  is  the 
dimension  of  NS{B).  | 


Since  we  have  a  matrix  operator  with  a  non-trivial  null  space,  it  is  conceivable 
that  in  solving  BS  =  /,  any  solution  o  may  have  components  in  NS{B).  As  we  are 
ultimately  reconstructing  an  image,  u  =  A'a,  the  effect  such  null  space  components 
have  on  the  appearance  of  the  reconstructed  image  is  of  great  importance.  We  will 
completely  characterize  NS{B)  by  constructing  a  basis  for  it. 

We  can  construct  a  basis  ,  9m-i]  for  NS(B)  in  the  following  fashion. 

All  of  the  qi  will  be  constant  by  angle,  with  the  constants  Qj  for  each  vector  defined 
as  follows: 
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\ 

ai 


02 


,  where  ai  G  TZ"' , 


^  j 


and 


ai 


(11  ...  If, 

(-1  -1  •••  -  If, 
(00  •••  Of, 


i  =  \ 
i  =  i  +  1 
otherwise 


Each  of  these  vectors  is  constant  by  angle,  and  Oj  =  0  for  each,  by  construction, 
so  qi  €  NS{B).  It  is  also  apparent  that  the  M-1  vectors  form  a  linearly  independent 
set.  Therefore  they  form  a  basis  for  NS (B), so  that  any  v  €  NS{B)  has  such  a 
representation.  Let  v  be  such  a  vector.  It  can  be  expressed  as 

M-l 

i=i 

This  expression  can  be  further  expanded  as 


\ 

(7i  +  72  +  •  •  •  +  7m-i)oi 


-710*2 


V 


-7203 


1 


— 7m-iOm 
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where  a,  G  7^"*  is  given  by 


a,  =  (1  1  if. 

It  is  clear  that  v  is  constant  by  angle,  and  that  these  constants  sum  to  zero. 

Any  vector  v  G  that  is  constant  by  angle,  (and  particularly  those  in 
NS(B))jdefines  the  coefficients  of  the  natural  pixel  representation  of  a  constant  image. 
Such  images,  when  displayed,  appear  as  uniform  shades  of  grey.  We  define  a  natural 
pixel  with  a  value  of  zero  to  be  black  ,  and  an  image  that  displays  to  uniform  black 
as  invisible. Beised  on  these  definitions,  we  have  the  following: 

Corollary  4.3:  If  v  ^  NS{B),  then  the  image  defined  by  u[x.,y)  =  A"v  is 
invisible. 

Proof:  Since  v  G  A^5(jB),  it  it  constant  by  angle,  and  the  sum  of  these 
constants  is  zero.  Now  u(a:,y)  =  A^v  =  Since  v  is  constant 

by  angle,  for  each  angle  i  the  natural  pixels  will  all  assume  the  same  value  to 
form  a  uniform  gray  sub-image  in  the  image  space.  Since  the  constants  that 
define  each  of  these  uniform  grey  sub-images  sums  to  zero,  the  combination 
of  the  sub-images  into  the  image  u(a:,  y)  will  be  a  uniform  image  in  the  image 
space  with  value  zero,  which  by  definition  is  invisible.  | 


These  results  are  significant.  Since  components  of  the  solution  in  the  NS{B) 
are  invisible,  we  need  not  concern  ourselves  with  them  in  the  framework  of  recon¬ 
structing  the  image.  Thus,  if  the  iterative  technique  chosen  to  solve  the  system  of 
equations  Ba  =  /  excites  components  in  NS(B)f^iS  far  as  the  display  of  the  image 
is  concerned  we  do  not  care.  (Such  components  do  not  affect  the  residual  calculation 
either). 

Our  goal  was  to  develop  a  discretization  that  resulted  in  a  linear  system  that 
could  be  solved  with  less  work  than  the  standard  square  pixel  discretization,  and 
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could  produce  comparable  results.  To  that  end,  we  now  compare  the  two  resulting 
linear  systems,  and  address  some  other  important  issues. 

C.  NATURAL  VS.  SQUARE  PIXELS 

In  general,  the  matrix  B  generated  by  using  natural  pixels  will  be  smaller  than 
that  generated  using  square  pixels,  for  the  same  source/detector  geometry.  If  we 
assume  M  angles  and  Ni{m)  detectors,  m  =  1  :  M,  then  N  =  Tlm=i  N\{m)  is  the 
number  of  rays  passing  through  the  image,  and  B  is  a,n  N  x  N  matrix.  Recall  that 
for  the  square  pixel  discretization,  we  divide  the  image  space  into  a  square  grid  of 
pixels.  The  resulting  matrix  K  is  of  size  N  x  n^.  Since  smaller  square  pixels  yield  a 
higher  resolution  image,  in  general  ri^  >  N  and  often  >>  N.  Thus,  the  matrix  B 
for  a  problem  discretized  by  natural  pixels  is  normally  smaller  than  the  matrix  K  for 
a  problem  discretized  by  square  pixels,  and  may  require  less  work  to  solve. 

On  the  other  hand,  K  is  very  sparse  w'hen  compared  to  B  .  Recall  that  non-zero 
entries  in  a  row  of  K  correspond  to  pixels  being  intersected  by  a  given  x-ray.  If  the 
width  of  the  ray  is  small,  say  on  the  order  of  the  width  of  a  pixel,  then  each  row  will 
contain  no  more  than  2n  non-zero  entries.  If  this  sparseness  is  exploited,  then  the 
work  needed  to  solve  the  system  can  be  greatly  reduced. 

Additionally,  K  is  in  general  a  rectangular  matrix,  and  the  only  practical  itera¬ 
tive  method  for  solving  the  resulting  system  of  equations  is  Kaczmarz’s  method  [Ref. 
17].  jB, however,  is  square  and  symmetric  and  possesses  all  of  the  other  properties 
discussed  above.  This  gives  us  more  selection  when  choosing  an  iterative  method  to 
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solve  the  accompanying  system  of  equations. 

One  very  significant  drawback  of  discretizing  by  natural  pixels  is  that  the 
natural  pixels  do  not  map  to  the  square  pixel  hardware  of  a  computer  screen  easily. 
One  way  to  display  the  image  is  now  outlined.  We  must  map  the  natural  pixels  to  a 
rectangular  grid  that  can  then  be  illuminated  on  a  computer  screen.  This  seems  at 
first  a  difficult  task,  but  it  can  be  easily  accomplished  through  the  use  of  a  square 
pixel  discretization  matrix  K  created  at  the  same  geometry  used  to  generate  B  . 

The  problem  we  are  attempting  to  solve  is  Ru  —  /,  where  R  is  the  Radon 
transform  operator.  Using  the  square  pixel  discretization,  we  discretize  both  the 
projection  and  image  spaces  and  arrive  at  a  linear  system  Kx  =  f.  When  we 
discretize  using  natural  |)ixels.  only  the  projection  space  is  discretized,  and  the  den¬ 
sity  function  u{x,y)  is  expanded  in  terms  of  strip  functions,  yielding  the  system 
Ba  =  AA'ol  —  f.  Nt)W.  if  we  view  x  as  an  approximation  of  A’ and  K  as  an 
approximation  for  A  ,  then  we  can  write 

X  =  K^a. 

Thus,  the  natural  pixels  can  be  mapped  to  square  pixels  just  by  backprojecting  a 
with  the  K  matrix. 

We  only  require  the  geometries  on  number  of  detectors  and  angles  to  agree 
for  both  discretizations.  The  natural  pixel  problem  is  independent  of  the  number  of 
square  pixels  to  which  it  is  being  mapped.  Therefore,  we  can  map  to  any  resolution, 
as  long  as  the  x-ray  geometries  agree. 
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The  cost  involved  to  display  the  image  is  primarly  due  to  building  the  matrix 
K  ,  which  can  be  excessive  if  high  resolution  is  required.  Fortunately,  this  is  a  one¬ 
time  cost,  for  K  can  be  repeatedly  used  to  display  different  images  generated  at  the 
same  geometry.  Another  drawback  is  that  K  might  be  too  large  to  fit  in  memory, 
so  storage  and  file-handling  procedures  will  be  required  when  displaying  the  image 
that  might  not  be  necessary  when  working  with  the  matrix  B  .  Also,  a  different  K 
is  required  for  each  geometry  used  to  generate  B  ,  or  for  different  resolutions  at  the 
same  geometry. 

We  have  developed  the  natural  pixel  discretization  for  the  image  reconstruc¬ 
tion  problem,  and  have  shown  that  it  produces  a  linear  system  whose  matrix  B  is 
symmetric  and  in  general  smaller  than  that  produced  by  the  conventional  square  pixel 
discretization.  An  analysis  of  this  matrix  produced  several  interesting  results.  The 
rank  of  B  is  a  function  of  the  x-ray  geometry  used  to  generate  it.  The  NS{B)  is 
characterized  by  vectors  that  are  constant  by  angle  with  constants  summing  to  zero, 
and  the  images  they  represent  are  invisible.  We  now  concentrate  on  solving  this  linear 
system. 
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V. 


GAUSS-SEIDEL  RELAXATION 


There  are  two  basic  strategies  for  solving  a  linear  system  of  equations  Ax  =  b, 
either  directly  or  iteratively.  The  direct  approach  generally  involves  a  factorization 
of  the  coefficient  matrix  A  .  Iterative  methods,  on  the  other  hand,  generate  a  se¬ 
quence  of  approximate  solutions  and  only  involve  the  matrix  A  for  matrix-vector 
multiplications.  If  A  is  large,  the  direct  approach  could  be  impractical  because  one 
direct  solve  costs  as  much  as  many  iteration  sweeps,  in  terms  of  floating  point  opera¬ 
tions.  The  effectiveness  of  an  iterative  method  is  determined  by  how  fast  the  sequence 
of  iterates  x^'^^  converges  [Ref.  19].  In  this  chapter  we  will  apply  the  Gauss-Seidel 
iterative  method  and  analyze  it  in  the  context  of  our  image  reconstruction  problem. 

A.  DERIVATION  AND  GENERAL  PROPERTIES 

We  derive  the  Gauss-Seidel  iterative  method  by  considering  how  to  solve  the 
linear  system  of  equations 

Qx  =  6, 

where  Q  is  a.n  N  x  N  matrix  assumed  to  be  nonsingular  with  nonzero  diagonal  entries. 
Then  the  formal  solution  is 

X  =  Q~^b. 
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Considering  a  splitting  of  the  matrix  Q  into  an  upper  triangular  part  U  ,  a  lower 
triangular  part  L  ,  and  a  diagonal  part  D  ,  so  that 

Q  =  D-L-U. 

We  can  now  write  the  system  of  equations  as 

Qx  =  {D-L-U)x  =  b, 


and  rearranging,  arrive  at 


{D  —  L)x  =  Ux  +  b. 


(5.1) 


We  convert  (5.1)  to  an  iterative  form  by  introducing  a  superscript  to  the  vector  x 


representing  its  index  in  the  sequence  of  approximations,  as  follows 


(/)- =  t/xW  +  6,  k>0. 


(5.2) 


Hence  the  {k  +  1)^'  approximation  is  generated  from  the  approximation,  and  the 
process  is  started  by  providing  an  initial  guess  Writing  (5.2)  in  terms  of  the 
elements  of  the  matrix  Q  produces 


f-i 


N 


j=l  j=t+l 


(5.3) 


which  can  be  further  rearranged  to  yield 


,.(‘+1)  _  —  I  6i  -  .  1  <  i  <  ^. 


(5.4) 


since  the  diagonal  elements  of  A  are  nonzero.  Equation  (5.4)  can  be  written  in 
algorithm  form  as 
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For  i  =  1,2,- N 


X 


(k+i) 


1 


j-i 


N 


\  j=\  7=H1  / 

which  defines  one  iteration  of  the  Gauss-Seidel  method.  Recasting  (5.4)  once  again 


in  matrix  notation,  we  get 


X 


=  (£)  -  +  (£>  -  L)-^6,  A:>0. 


(5.5) 


Defining  Pq  =  {D  —  L)  and  c  =  {D  —  L)  *5,  the  Gauss-Seidel  method  (5.5) 
becomes 


=  Pgx^’^^  +  c,  A->0. 


The  matrix  Pc  is  known  as  the  Gauss-St  idd  ilt  ration  matrix  . 

Convergence  of  the  Gauss-Seidel  nH‘lho<l  ran  be  analyzed  in  IcTins  of  the  error 
err/or, defined  as 


where  x*  is  the  exact  solution  of  the  proldi'in.  Sinrt*  /*  i.s  the  exact  solution,  then 
X*  =  Pqx*  +  c,  so  that 

=  n;/*'4r-(n;r +  r) 

=  Pa(P^'^-r) 

=  Pgc^^K 

Thus  we  arrive  at  the  relationship 

g{k+i)  ^  (5.6) 
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Recursive  applications  of  (5.6)  yields 


=  •••  = 

Thus,  we  get  convergence  if  0  as  fc  — >  oo,  and  this  convergence  is  inde¬ 

pendent  of  the  initial  guess 

Define  the  spectral  radius  of  Pq  =  p{Pg)  to  be  the  magnitude  of  the  largest 
eigenvalue  of  Pq.  The  following  well-known  theorems  [Ref.  22]  give  conditions  under 
which  the  Gauss-Seidel  method  is  guaranteed  to  converge: 

Theorem  5.1:  (Pg)*^^^  0  as  k  —y  oo  if  and  only  if  p{Pg)  <  1- 

Theorem  5.2:  If  A  £  symmetric  and  positive  definite,  and  Pa  is 

the  Gauss-Seidel  iteration  matrix  formed  from  A,  then  p{Pg)  <  1- 

Thus,  a  condition  that  ensures  convergence  of  the  Gauss-Seidel  method  to 
is  that  the  matrix  defining  the  system  of  equations  to  be  solved  be  symmetric  and 
positive  definite.  The  matrix  B  from  the  natural  pixel  discretization  is  symmetric  and 
positive  semi-definite,  so,  apparently,  Gauss-Seidel  applied  to  it  is  not  guaranteed  to 
converge.  However,  it  will  be  shown  that  Gauss-Seidel  is  convergent  for  the  matrix 
P,  and  why  this  must  be  so. 

B.  CONVERGENCE  ANALYSIS 

First,  we  show  that  the  Gauss-Seidel  method  cannot  diverge  for  this  problem. 
We  introduce  the  following  definitions: 
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Definition  5.1:  The  energy  inner  product  can  be  defined  as  (x,y)  =  {Gx,y), 
where  G  is  semidefinite. 


Definition  5.2:  The  energy  norm  is  defined  as  |||x|||  = 


The  energy  norm  is  actually  a  semi-norm  ,  because  |||x|||  =  0  does  not  imply 
that  X  =  0.  We  will  show  that  when  measured  with  the  energy  norm,  Gauss-Seidel 
cannot  diverge. 

Recall  that  (5.4)  is 


X 


(k+i)  _  1 


t-1 


=  -b.-E'j.-r"-  E 


N 


.(fc) 


I  <i<  N, 


i=i 


j=t+i 


and  after  further  manipulation  we  obtain 


.(fc+i) 


1-1 


N 


.w 


Qi 


(5.7) 


J=1 


Let 


j=i  j=i+i 

be  the  component  of  the  current  residual  vector  f*  during  the  step  of  the 
iteration.  That  is,  it  is  computed  using  both  and  and  changes  for  each 

value  of  i.  Substituting  into  (5.7)  results  in 

xj'^+i)  =  3.W  _  II,  1  <  z'  <  yv.  (5.8) 

qa 

Let  Wi  to  be  the  standard  basis  vector  for  that  is,  the  vector  with  a  1  in 

the  position  and  zeros  in  all  other  positions.  Then  we  can  express  the  diagonal 
elements  of  Q  as 

qu  =  wi'^Qwi  =  {Qwi,wi) 
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and  r,  as 


n  =  f^Wi  =  {f,Wi). 


Writing  (5.8)  in  terms  of  these  inner  products  produces 


•*'1 


{r,Wi) 

{Qwi,wi)' 


I  <i  <  N. 


(5.9) 


Only  one  component,  the  of  the  approximation  is  being  changed  at  each  step,  so 
that  (5.9)  can  be  written  in  terms  of  vectors  cis 


f(fc+i)  ^  w,.  1  <  /  <  .V.  (5.10) 

{Qw,,w,) 

where  the  arrow  indicates  replacement,  or  overwriting.  Tin*  above  formulations  are 
used  to  yield  the  following  result: 


Theorem  5.3:  The  Gauss-Seidd  tiitlhod  opplnd  to  lin  =  /.  u'hen  B  is  the 
matrix  generated  via  discretization  bg  natural  pij(l.'>.  ».*•  boundtd  in  the  energy 
norm. 

Proof:  Define  the  error  vector  as 


(5.11) 

where  <3*  is  an  exact  solution,  at>d  m)te  tliat 

=  B(Q’-a(''))  =  /y<r  - =  f-Bo^^^  =  f.  (5.12) 
Taking  the  energy  norm  of  the  error  vector  gives 


r(fc+i)|||2  _ 


..|e' . Ill-  = 

and  applying  (5.10),  (5.11)  and  (5.12)  to  the  right  side  produces 


'gW  _  <  ggf  .f.  > 


<  Wi  > 
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Some  algebraic  manipulations  produce  the  following 

\  V  <Bwi,Wi>  , 


-Bwi,  \e'‘  ’ - — i: - Wi 


<  Bwi,Wi  > 


<  Bwi,wi  > 


which  can  be  expanded  as 


|e'‘+"||P  =  <  Be-W.e-W  >  ^  <  Bel‘’\w,  > 

<  Bwi,Wi  > 

<Be^^\wi>  <Be^^\wi>'^  ^ 

<  Bwi,e^^’  >  H - - —  <  Bwi,Wi  > 


<  Bwi,  W{  > 
Simplification  results  in 


<  Bwi,Wi  >2 


=  <  BeW.eW  >  -2< 

<  Bwi^Wi  >  <  Bwi,Wi  > 

=  <  B?W,e-W  >  -  < 

<  Bwi,  Wi  > 


=  e 


(fc)|||2  _ 


<  Be^^\  Wt 


<  Bwi,  Wi  > 

Since  B  is  positive  semidefinite  and 

<  Bwt^Wi  >  =  bii  >  0, 

we  have 

|||J(W)|||2  <  |||e-'‘)||p. 
Therefore,  the  norms  remain  bounded. 


This  useful  result  assures  us  that  the  Gauss-Seidel  iteration  can  be  applied  to 
Ba  =  /,  and  provided  we  use  the  energy  norm  to  measure  performance  we  need 
not  fear  a  divergent  process.  To  obtain  a  stronger  (convergent)  result,  we  show  that 
under  certain  conditions  Gauss-Seidel  is  equivalent  to  Kaczmarz’s  method,  and  that 
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the  convergence  of  Kaczmarz  established  in  Theorem  3.1  carries  over  to  Gauss-Seidel 
applied  to  the  natural  pixel  problem. 

Consider  applying  Kaczmarz’s  method  to  Au  =  /,  where  A  :  H  , 

u  £  H ,  f  £  71^ ,  H  is  a.  real  Hilbert  space  spanned  by  the  constant-valued  strip 
functions,  and  A  is  defined  (see  page  47)  so  that 


\ 


Au 


{rp2,u) 


If  €  Range(/1*),  then  the  sequence  generated  by  Kaczmarz’s  method  con¬ 
verges  to  the  minimum  L2  norm  solution  as  k  00.  One  sweep  of  Kaczmarz  can  be 
expressed  as 

Set  u  = 

For  i  =  1,2, • • • , N 

Solve  <  ty,,  A{u  -f-  sA*Wi)  —  f  >=  0  for  s 
Set  u  =  u  -f  sA'*Wi. 

Set  =  u 

Solving  for  s  in  the  above  algorithm 

0  =  <  Wi,A{u  + sA*Wi)  —  f  > 

=  <  Wi,  Au  -f-  sAA*Wi  —  f  > 
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=  Wi^  Au  +  swi^  AA*Wi  —  Wi^  f 

=  <U,  ijji  >  +5  <  V’t-  >  -fi, 


yielding 

^  _  fi  -  {i’i-.u) 

<  ■0,-,  0t-  >  * 

Thus,  Kaczmarz’s  method  becomes 


Set  u  = 


For  i  =  1, 2, • • • , 


u  =  u  -\- 


fi  -  ii^uu) 


Set  =  u 


Now,  consider  applying  the  Gauss-Seidel  method  to  the  natural  pixel  problem 
lie!  =  /,  where  B  =  A  A'  and  u  =  .4*0.  Note  that  u  €  Range(.4’).  Gauss-Seidel 
1)11  tills  problem  can  be  written  as 
Set  Q  = 

For  i  =■  1,2,  •  •  • , 

Solve  <  in,-,  B[a  +  sin,)  —  /  >=  0  for  s 
Set  <3  =  a  +  swi. 

Set  =  a. 

Again,  we  can  solve  for  s  in  the  above  algorithm.  Let 
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B  = 


-T  \ 

h. 


-T 


,  where  6,-  =  pa  ■  ■  ■  Atv). 


Then,  solving  for  s  ,  we  obtain 


yielding 


But 


0  =  <  Wi,  B(a  +  swi)  —  f  > 

=  <  Wi,  Ba  +  sBwi  —  f  > 

=  Wi^  B  a  +  swi^  Bwi  —  Wi^  f 

-T 

=  b,  a  +  spa  -  fi, 


hi  a  =  <  A‘a,xl>t  >  =  <  u,xl)i>  and  ;9„  =  <  >? 


so 


_  fi-  <  U,xJ)i  > 

<tl)i,'4)i>  ' 

and  after  substituting  into  the  Gauss-Seidel  algorithm  we  can  write 


<3  =  <3  + 


fi-  <U,xf)i>  ^ 

- ; — - Wi- 

<  tPi,  Vi  > 
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Multipliying  both  sides  by  A'*  yields 


A  o.  =  A  a  -\ - ; — - . 


But  A*o.  =  u,  so  we  arrive  at 


u  =  u  + 


fi  -  (V’P  u)  I 

{i’i,  i’i) 


Therefore  Gauss-Seidel  applied  to  Ba  =  /  is  equivalent  to  Kaczinarz  vi^plied 
to  Au  =  /  [Ref.  2].  This  observation  yields  the  following  result: 


Theorem  5.4:  If  £  Range(B^ ),  then  the  sequence  generated  by  the 
Gauss-Seidel  method  converges  to  a  such  that  it  =  A*6(  is  the  minimum  J  -i 
norm  solution  of  Au  =  f  as  k  oo,  provided  such  a  solution  exists. 

Proof:  Assume  Ba  =  f  is  consistent.  If  €  Range(B^),  then 

BaW  =  AA’a^^K 

so  we  can  define  £  Range(A*).  By  Theorem  3.1,  we  know  that 

the  sequence  generated  by  Kaczmarz’s  method  converges  to  the  minimum 
Z/2  norm  solution,  u,  of  Au  =  f  as  k  —i’  oo.  But  Kaczmarz  is  equivalent  i  > 
Gauss-Seidel  under  these  conditions,  so  Gauss-Seidel  must  converge  to  some  o 
such  that  u  =  A'a.  | 


Next,  we  show  that  the  problem  Au  =  /  can  always  be  made  to  be  consistent 
so  that  a  solution  exists,  and  that  we  can  always  produce  an  inital  guess  in  the 
Range{B^).  Consider  a  projection  of  the  image  at  a  given  angle,  that  is,  the  collection 
of  emergent  x-ray  intensities  measured  by  the  array  of  detectors.  The  natural  pixel 
discretization  assumes  complete  coverage  of  the  image  by  x-rays  at  each  angle.  Noting 
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that  the  density  of  the  image  remains  unchanged  from  angle  to  angle,  and  that  the 
initial  x-ray  intensities  are  constant  from  angle  to  angle,  then  the  total  emergent  x-ray 
intensities  for  each  projection  must  be  constant  as  well.  This  must  be  the  case  as  the 
image  will  absorb  the  same  amount  of  energy  independent  of  the  direction  of  that 
energy.  We  say  that  the  data  /  is  compatible  with  the  image  reconstruction  problem 
when  the  total  emergent  x-ray  intensities  are  constant  for  each  angle,  that  is,  the  sum 
of  the  elements  in  the  piece  of  /  corresponding  to  each  angle,  is  constant.  Given  this 
definition,  we  have  the  following 

Theroem  5.5:  If  the  data  f  is  compatible  with  the  image  reconstruction  prob¬ 
lem,  then  Au  =  /  is  consistent  and  has  a  solution. 

Proof:  Let  /  be  compatible.  Then  /  ^  NS{B)  because  it  cannot  be  constant 
by  angle  with  the  Qi  =  0  and  compatible  at  the  same  time.  Therefore 
/  ^  NS{A“').  Assume  the  sum  of  the  elements  in  /  corresponding  to  each  angle 
sum  to  the  constant  C.  Let  q,  be  the  i*^  basis  vector  for  NS(B)  as  constructed 
in  Chapter  IV,  and  consider  <  f,qi  >.  If  the  inner  product  is  computed  by 
angle,  the  resulting  sum  is 

^  I  iQt  ^  —  G-fO  +  O-l-  •••  d-O  —  C-fO-|-  •••  -fO  =  0, 

where  the  negative  entry  is  in  the  (i  -f  1)*‘  position.  Therefore  /  is  orthogonal 
to  each  basis  vector  of  NS(B)  and  hence  f  J.  NS{B)  =  NS{A').  Thus 

/  €  Range{A).  | 


If  for  some  rezison  the  data  we  are  given  to  reconstruct  is  not  compatible  due 
to  measurement  errors,  we  can  correct  the  data  and  enforce  compatibility  by  adding 
some  constant  Ci  to  each  of  the  elements  of  the  view.  Thus  we  can  always  ensure 
a  consistent  system,  and  theoretically  Gauss-Seidel  will  converge. 

Next,  we  use  the  concept  of  compatibilty  to  define  an  initial  guess  for  the 
Gauss-Seidel  iteration.  Recall  that  if  /  is  compatible,  then  the  sum  of  its  elements 
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over  any  angle  is  constant,  call  it  C.  Consider  a  vector  that  is  constant  by  angle, 
where  the  constants  are  such  that 

M 

c  =  x;®.- 

t=l 

We  define  such  a  vector  to  be  the  average  grey  representation  of  the  image  whose 
data  vector  is  /.  There  is  an  infinite  set  of  constants  a,-  that  can  define  an  average 
grey  representation  of  a  given  image.  We  show  that  at  least  one  of  them  defines  a 
vector  in  the  Range{B^). 


Lemma  5.1:  Let  v  €  be  an  average  grey  rrpn s(  niation  of  some  image 
u(Xfy)  with  constants  a,-  such  that 


n. 


for  all  i.  Then  v  G  Range{B^). 


£ 

M 


Proof:  Let  a  be  the  arithmetic  mean  of  the  ron>taTi!s  o,  in  llir  definition  of 
constant  by  angle.  Now  v  can  b(‘  written  a> 


V  =  t  H  :  »  \ . 

where  Vfi  G  Range{B^)  and  €  \S[B).  \\r  claim  this  decomposition  can 
be  expressed  as 


o’,  \ 

n  \ 

f  at  -  it  y 

-• 

OL2 

Cl 

Cl;  —  Cl 

. 

- 

oTm  ) 

ith  -  it  ) 

where  the  quantities  a*  are  constant  vectors  of  length  !^\{i)  with  value  a,.  A 
well  known  theorem  from  statistics  [Ref.  2*1]  states  that  i)  =  0. 

Therefore,  by  Theorem  4.3  the  vector 


VN  = 


(  cTi  —  a  \ 

\  a  M  -  ^  J 
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is  constant  by  angle  and  in  NS(B).  But  we  have 

Em  ^ 

_  i=l  M  ^  il 

MM' 

Therefore  iT;v  =  0  and  v  =  vr.  Thus  v  G  Range{B^).  | 

We  have  shown  that  the  Gauss-Seidel  method  applied  to  the  natural  pixel 
problem  is  equivalent  to  the  method  of  Kaczmarz  applied  to  the  square  pixel  problem, 
and  given  an  appropriate  initial  guess  Gauss-Seidel  converges  to  the  same  minimum 
Z/2  norm  solution  as  Kaczmarz.  Furthermore,  we  have  shown  roiiditioiis  such  that  the 
problem  is  consistent  and  determined  specifically  what  an  ajipropriate  initial  guess 
should  be.  Analysis  of  the  spectral  properties  of  the  Gauss-Seidel  iteration  matrix 
Pg  is  still  required,  to  gain  further  insight  as  to  why  the  method  converges,  and  to 
analyze  the  rate  of  convergence. 

C.  SPECTRAL  ANALYSIS 

It  hcis  been  shown  thus  far  that  the  Ganss-Seiilel  method,  when  applied  to 
the  linear  system  derived  from  the  natural  pixel  discretization,  remains  bounded 
ill  the  energy  semi-norm  and  that  for  rertain  starting  vectors  it  converges  to  the 
minimum  L2  norm  solution.  The  matrix  H  is  only  positive  semi-definite,  so  Theorem 
5.2  does  not  apply  and  the  Gauss-Seidel  iteration  matrix  Pc  is  not  guaranteed  to 
have  a  spectral  radius  p{Pg)  <  1-  Non-divergence  has  been  established  by  Theorem 
5.3,  implying  p{Pg)  ^1-  If  an  appropriate  starting  vector  is  given.  Theorem  5.4 
guarantees  convergence,  implying  that  the  eigenvectors  corresponding  to  eigenvalues 
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A,-  of  Pa  with  |A,|  =  1  do  not  affect  the  norm  of  the  residuals.  If  such  components 
affected  the  measure  of  the  iterates,  then  the  error  expression 


would  not  decrease  with  k  .  We  use  this  fact  to  arrive  at  the  following  result: 

Theorem  5.6:  If  zl  is  an  eigenvector  of  Pq,  and  if  Zi  6  NS(B),  then  the 
corresponding  eigenvalue  A,-  =  1. 

Proof:  For  the  Gauss-Seidel  method  applied  to  the  linear  system 

B5  =  /, 

we  have 

B  =  D  -U  -  L,  Pg  =  {D-L)-^U, 

and 

Let  an  exact  solution  to  the  system  be  o'.  Then  the  iteration  becomes 

=  {D-  +  {D-  L)-^Ba. 

Writing  this  expression  in  ternas  of  the  error  vector  yields 

=  PgC^^^ 

=  (D  - 

=  {D  -  L)-\D  -  -  {D  -  L)-^  Be^^^ 

=  -  (D  - 

Now,  let  e  =  Zj  G  NS{B).  Then 

=  Zi-{D  -  L)-^Bzi  =  Zi. 

But  we  know  that  =  Pae^^'^  as  well,  so 

g(fc+i)  ^  p^^(k)  ^  ^  XiZi  =  Zi  =  e^'^\ 

Therefore,  A,  =  1.  | 
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We  have  learned  that  even  though  the  spectral  radius  of  the  Gauss-Seidel 
iteration  matrix  for  this  problem  may  equal  unity,  the  method  still  converges  because 
the  eigenvectors  of  Pq  corresponding  to  eigenvalues  of  modulus  one  are  in  NS(B). 
Therefore,  they  do  not  contribute  to  the  approximation.  Theoretically,  Gauss-Seidel 
should  be  able  to  solve  the  natural  pixel  problem.  What  remains  to  be  done  is  actually 
apply  Gauss-Seidel  to  several  linear  systems  and  analyze  the  resulting  behavior. 

D.  BEHAVIOR  OF  GAUSS-SEIDEL  APPLIED  TO  THE 
PROBLEM 

Gauss-Seidel  is  applied  to  several  linear  systems  generated  at  varying  geome¬ 
tries,  and  with  right-hand-side  data  created  both  analytically,  and  by  projecting 
computer-generated  images.  In  all  cases,  the  overall  behavior  observed  is  of  rapid 
initial  convergence  that  eventually  stalls  out,  just  as  occurs  with  the  Kaczmarz  iter¬ 
ation.  Figure  33  depicts  plots  of  both  the  norm  of  the  residual  and  the  convergence 
factor  plotted  against  iterations  for  a  typical  problem.  Here,  the  convergence  factor  is 
formed  as  the  ratio  of  successive  residual  norms.  The  magnitude  of  the  residual  error 
in  all  cases  is  on  the  order  of  10~^  or  10“^,  well  short  of  machine  precision.  To  explain 
why  this  behavior  occurs,  the  spectrum  of  the  matrix  B  is  analyzed.  Even  though  B 
is  square,  to  parallel  the  analysis  presented  in  Chapter  III  we  use  the  singular  value 
decomposition. 

Recall  the  singular  value  decomposition  (SVD)  of  B 

B  =  (/SV^, 
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Norm  of  Rouduai  vs  Itofabons  Convergence  Factor  vs,  Iteralons 


Figure  33.  Convergence  of  Gauss-Seidel  on  a  typical  problem. 


where  U  and  V  are  orthogonal  matrices,  and  E  is  a  diagonal  matrix  whose  diagonal 
entries  a,  are  the  singular  values  of  B  .  The  columns  of  U  and  V  are  known  as  left 
and  right  singular  vectors  ,  respectively.  Singular  values  are  real,  nonnegative,  and 
ordered  such  that 

•  >  <7r  >  0  •  •  •  0. 

The  number  of  non-zero  singular  values  r  <  N  equals  the  rank  of  B  .  The  SVD  can 
be  rewritten  eis 

BV  =  t/E, 

and  if  the  columns  of  this  expression  are  compared,  we  arrive  at  a  collection  of  linear 
systems 

Bvi  =  aiU{,  1  <  i  <  N, 

where  Ui  and  Vi  denote  the  columns  of  U  and  V,  respectively.  It  is  possible  to 
determine  which  singular  values  have  singular  vectors  Vj  that  are  slow  to  converge. 
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Figure  34.  Singular  Values  for  32  Detectors  at  8  and  20  Angles. 


This  is  accomplished  by  using  Gauss-Seidel  to  solve,  for  each  z,  the  problem  Ba  — 
OiUi.,  using  zero  as  the  initial  guess.  The  columns  of  V  form  an  orthonormal  basis  for 
.  so  any  image  can  be  expressed  as  a  linear  combination  of  the  If  the  image 
has  components  which  are  slow  to  converge,  then  the  iteration  will  stall.  Figure  34 
depicts  the  spectrum  of  singular  values  for  two  geometries.  One  is  constructed  using 
32  detectors  at  each  of  8  evenly  spaced  angles,  producing  a  matrix  of  size  232  x  232 
and  of  rank  225.  The  second  geometry  is  32  detectors  at  each  of  20  evenly  spaced 
angles,  resulting  in  a  matrix  of  size  592  x  592  and  of  rank  573.  It  should  be  pointed 
out  that  for  a  geometry  of  20  angles  and  32  detectors  per  angle,  one  would  think 
there  would  be  640  total  x-rays  producing  a  matrix  B  of  size  640  x  640.  However, 
only  592  x-rays  actually  pass  through  the  image  square,  which  explains  the  disparity 
in  the  size  of  the  matrix.  This  fact  holds  for  all  geometries. 

Notice  that  in  both  graphs,  the  magnitude  of  the  singular  values  can  be  divided 
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into  three  bands.The  first  of  these  is  the  left  portion  of  the  plots  where  the  curve  is 
nearly  horizontal,  the  second  is  in  the  center  where  the  magnitude  of  the  singular 
values  noticeably  decrease,  and  the  third  is  to  the  right  where  the  singular  values 
drop  suddenly  to  zero  -  the  null  space.  As  before,  we  will  define  the  center  band  as 
the  near  null  space^dind  the  left  band  as  the  resolvable  regionJt  will  be  observed  that 
singular  vectors  corresponding  to  singular  values  in  the  near  null  space  are  the  slow 
components  of  the  image  to  converge. 

Consider  using  Gauss-Seidel  to  solve  the  SVD  system 

BVi  =  (TiUi, 

whose  exact  solution  is  tT,,  for  various  values  of  i .  The  resulting  solutions  Vi  should  be 
good  apj)roximations  for  the  corresponding  uj.  We  can  decompose  the  Vi  into  linear 

combinations  of  the  singular  vector  basis  as 

N 

Vi  =  Y^PjVj. 

J=i 

The  singular  vector  basis  is  orthonormal,  and  the  0j's  can  be  computed  as 

/3j  =  vi^vj,  1  <  i  <  A^. 

Hence  we  can  determine  the  components  of  Vi  in  the  directions  of  each  of  the  singular 
vectors  Vj.  For  the  exact  solution  uj,  we  have 

1,  i  =  j 

0,  otherwise. 
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Spoctral  Docompositon  (or  Singular  Valuo  #175  Spoctral  Dacompoeiton  for  Singular  VaJu«  #525 


Figure  35.  Approximations  for  Singular  Values  in  the  Resolvable  Region  and  Near 
Null  Space. 

A  plot  of  the  (dj  for  the  exact  solution  will  be  a  spike  of  magnitude  one  at  index  r.The 
following  figures  are  plots  of  the  \(3j\  for  some  indices  both  in  the  resolvable  region 
and  in  the  near  null  spare  for  a  geometry  of  32  detectors  at  20  angles. 

A  qualitative  interpretation  of  these  plots  follows.  Three  things  are  readily 
apparent.  First,  components  in  the  resolvable  region  are  almost  totally  recovered 
(Figure  35,  left).  A  spike  of  magnitude  one  located  at  the  appropriate  index  is  clearly 
present,  along  with  a  small  amount  of  noise.  These  components  do  not  adversely  affect 
the  performance  of  the  iteration.  Next,  components  in  the  near  null  space  are  not 
recovered  well  (Figure  35,  right).  There  is  a  partially  recovered  spike,  and  significant 
noise  is  present.  These  components  represent  the  unrecoverable, or  s/ow, components 
of  the  image  that  cause  the  iteration  to  stall.  Finally,  iterating  with  Gauss-Seidel, 
just  as  with  Kaczmarz,  mixes  modes, that  is,  it  introduces  additional  components  of 
the  singular  vector  spectrum  as  noise  into  the  approximation  vi.  The  noise  includes 
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components  in  the  null  space,  but  is  predominantly  composed  of  near  null  space 
components,  that  are  not  part  of  the  actual  solution 

This  last  observation  is  significant  in  that  we  can  have  an  image  that  lives 
entirely  in  the  resolvable  region,  apply  Gauss-Seidel  to  reconstruct  it,  and  excite 
components  in  the  approximation  that  live  in  near  null  space  and  null  space.  The 
components  in  the  null  space  are  not  a  problem,  as  the  images  they  generate  have 
been  shown  to  be  invisible,  and  they  do  not  affect  the  measurement  of  the  residual. 
However,  those  in  the  near  null  space  are  a  problem,  for  as  seen  above,  they  are  slow 
to  be  recovered  and  virtually  all  of  the  residual  error  can  be  attributed  to  them. 

In  spite  of  these  apparent  shortcomings,  Gauss-Seidel  applied  to  the  natural 
pixel  discretized  problem  reconstructs  images  quite  well.  The  following  figures  depict 
actual  and  reconstructed  images  for  a  brain  phantom,  ceated  by  superimposing  a 
collection  of  ellipses  and  rectangles  of  varying  grey  levels  on  each  other.  The  data 
vector  /  W21S  then  generated  by  projecting  the  image  with  Kaczmarz  matrices  of 
assorted  geometries. 

The  behavior  described  above  is  not  unique  to  the  Gauss-Seidel  method.  Under 
the  assumption  that  other  types  of  iterative  methods  might  not  exhibit  the  behavior 
of  Gauss-Seidel,  several  Lanczos-based  methods  were  applied  to  the  problem  as  well. 
In  particular,  the  algorithms  SYMMLQ  and  MINRES  [Ref.  25]  and  several  incom¬ 
plete  orthogonalization  methods  [Ref.  26]  were  applied  to  the  image  reconstruction 
problem.  In  all  cases,  convergence  was  initially  rapid,  followed  by  stalling.  The  recon- 
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Figure  36.  Actual  and  Reconstructed  Images  -  Brain  Phantom 


structed  images  produced  by  these  methods  were  no  better  than  those  produced  by 
Gauss-Seidel.  As  all  of  these  methods  require  more  work  per  sweep  than  Gauss-Seidel 
and  produce  no  better  results,  we  will  not  examine  them  further. 

We  have  conducted  an  analysis  of  the  Gauss-Seidel  method  as  applied  to  our 
image  reconstruction  problem.  Although  convergence  ha^  been  established,  there 
are  certain  components  of  the  solution  that  are  slow  to  be  recovered.  In  the  next 
chapter,  a  multilevel  method  will  be  developed  that  accelerates  the  convergence  of 
the  Gauss-Seidel  iteration. 
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VI. 


A  MULTILEVEL  APPROACH 


A.  BASIC  MULTILEVEL  CONCEPTS 

Multilevel  methods  were  developed  to  overcome  the  numerical  stalling  of  iter¬ 
ative  methods  for  numerical  partial  differential  equations.  A  fundamental  principle 
of  any  multilevel  methodology  is  that  the  amount  of  computational  work  should  be 
proi)ortional  to  the  amount  of  real  physical  changes  in  the  computed  system  [Ref.  27]. 
If  it  is  not,  e.g.  if  successive  sweeps  of  an  iterative  method  on  a  linear  system  produce 
.smaller  and  smaller  reductions  in  the  error,  then  a  more  efficient  method  should  exist 
to  ap|)roach  the  problem.  The  image  reconstruction  problem  exhibits  such  behavior. 
I'lie  situation  occurs  when  there  exists  several  solution  components  with  different 
scales  that  conflict  with  each  other.  The  answer  could  be  a  multilevel  approach, 
which  involves  interactively  employing  several  scales  of  discretization  to  resolve  such 
conflicts,  avoid  stalling  of  the  iteration,  and  elintinate  computational  waiste.  Before 
such  a  method  can  be  designed,  it  is  necessary  to  understand  the  basic  methodology 
behind  the  approach  and  why  it  works.  To  that  end,  a  coarse  grid  correction  scheme 
as  applied  to  a  simple  one-dimensional  partial  differential  equation  will  be  developed 
and  analyzed. 
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1.  Elements  of  a  Multilevel  Method 

Consider  the  second-order  differential  equation  for  the  steady-state  heat 
distribution  in  a  rod  of  uniform  density  that  has  the  temperature  fixed  at  both  ends 

—u"{x)  —  au{x)  =  f{x),  0<a:<l,  cr>0 

subject  to  the  Dirichlet  boundary  conditions 

u(0)  =  u(l)  =  0. 


The  problem  is  discretized  by  breaking  its  domain  into  equal  subin¬ 
tervals  of  width  h  ,  which  defines  the  node  points  Xj  =  jli.  j  =  0,  !,•••,  A'^, 

and  forms  a  grid  which  will  be  denoted  as  iV‘.  If  wt*  let  approximate  u{xj),  and 
if  we  approximate  u''{xj)  with  a  finite  2’“^  order  «li(rerenre.  then  tlie  problem  can  be 
rewritten  as 


-Uj_i  d-  2vj  -  Vj+i 

h? 


+  avj  =  J{x,). 


=  r\  =  U.  j  =  I  :  A'  —  1 . 


Here,  we  have  converted  an  ordinary  <liff«T«-ntial  (‘<|uatiun  into  a  system  of  N-1  alge¬ 
braic  equations  in  N-1  unknowns,  with  the  err«»r  in  the  approximation  being  of  0{h}). 
If  we  define 


X  = 


>1 

(  \ 

(  \ 

\ 

Xi 

V\ 

/U.) 

u{xi) 

X2 

V2 

U{X2) 

; 

,  V  = 

: 

,  /  = 

,  and  u  = 

: 

a:jV-i  j 

^  VN-1  j 

^  f{xN-l)  y 

^  u{xN-i)  ^ 
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then  we  can  express  the  problem  as  the  linear  system 


where 


2  + /i^cr  —1 

—  1  2  +  h'^a  — 1 


\ 


9 


-1 


—  1  2  +  h^a  j 


and  whose  solution  v  approximates  u  =  u(x)  at  the  grid  points.  It  can  be  determined 
analytically  that  the  eigenvectors  of  A  are 


-t-  = 


sra(^) 


sin(2^) 


,  A-  =  I  :  A'  - 


Graphing  one  finds  that  the  graph  is  smooth  for  small  values  of  k  ,  and  becomes 
increasingly  oscillatory  as  k  increases,  as  shown  if  Figure  37.  Following  the  analysis  of 
the  previous  chapter,  we  consider  applying  an  iterative  method  to  this  linear  system. 
We  will  use  the  weighted  Jacobi  method  as  a  vehicle  for  this  discussion.  Letting 


A  =  D-L-U, 
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then  the  weighted  Jacobi  iteration  matrix  is  given  by 


Pj  =  {l-oj)I  +  uD-^{L  +  U), 

where  u  G  71  is  a.  weighting  factor  to  be  chosen.  It  can  be  shown  that  the  eigenvectors 
of  the  matrix  Pj  are  identical  to  the  eigenvectors  of  the  matrix  A  .  The  eigenvalues 
of  Pj  are 

aTT 

At(Fj)  =  l-2wsin"(— ),  fc  = 


U1 

Figure  37.  A  mode  at  k=^  ou  grids  of  N=12  and  N=6. 
We  define  the  error  vector  as 


e  =  V  —  u, 


(6.1) 


so  that  the  initial  error  in  the  iteration  is  e  This  can  be  expanded  in  terms  of  the 
eigenvectors  of  Pj  as 

N-l 

g'(o)  _  ^ 

k=i 
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After  n  sweeps  of  the  iteration,  we  have 


=  E  =  II  <^kKiPj>'’k. 

k=l  k=\ 

Thus,  after  n  sweeps  the  error  has  been  reduced  by  a  factor  of  A^(Pj).  This  means 
that  for  large  values  of  k  ,  the  error  is  reduced  rapidly,  while  for  small  values  of  k 
it  is  not.  We  will  call  the  eigenvectors  with  index  1  <  A;  <  y  the  low-frequency  or 
smooth  modes,  and  those  with  index  y  <  A:  <  —  1  the  high-frequency  or  oscillatory 

modes.  Our  analysis  shows  that  error  components  corresponding  to  smooth  modes 
are  slow  to  be  eliminated.  This  behavior,  while  not  as  easily  quantified,  holds  for  the 
Gauss-Seidel  method  as  well  (Ref.  28].  The  aim  of  a  multilevel  approach  is  to  devise 
a  way  to  address  these  slow  components. 

To  develoj)  the  coarse  grid  correction  scheme,  we  first  must  define  a 
second  set  of  y  grid  points  obtained  by  selecting  every  other  grid  point  from 
17^.  Note  that  17^^  is  coarser  that  meaning  that  the  grid  spacing  is  wider.  Now, 
cLSsume  the  iterative  method  has  been  applied  to  the  linear  system  on  the  original 
grid  17^,  until  only  smooth  error  components  remain.  If  we  consider  what  the  smooth 
error  components  look  like  on  17^^,  as  illustrated  in  the  following  Figure  38,  we  see 
that  on  the  coarser  grid  the  error  becomes  more  oscillatory.  On  f7^  mode  4  is  one 
third  of  the  way  up  the  spectrum.  However,  on  i7^^  mode  4  is  now  two  thirds  of  the 
way  up  the  spectrum  and,  therefore,  more  oscillatory. 

To  be  more  precise,  consider  the  mode  on  f7^  evaluated  at  the  even- 
numbered  grid  points,  which  is  exactly  f7^^.  We  can  write 
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(6.2) 


^  — 
Zu  - 


2/c7r  ^ 


sm(^) 


sin(^) 


sin(4«^) 


where  the  superscripts  denote  the  two  grids.  This  implies  that  the  mode  on 
becomes  the  mode  on  as  long  as  1  <  A:  <  Thus  in  moving  from  the  finer 
grid  to  the  coarser  grid,  a  mode  becomes  more  oscillatory,  and  as  such  relaxation 
should  be  more  effective.  It  is  equally  important  to  note  that  the  smoothness  of  e 
after  relaxation  on  is  what  allows  us  to  go  to  -  only  if  e  is  smooth  can  it  be 
accurately  represented  on  a  coarser  grid. 


Figure  38.  j4  mode  at  k=4  grids  of  N=12  and  N=6. 

It  can  also  be  shown  that  if  A;  =  ^ ,  the  mode  becomes  the  zero  vector 
on  and  ii  y  —  ^  ^  ~  mode  on  becomes  the  (N  —  k)^^  mode 

on  This  last  statement  says  that  an  oscillatory  mode  on  the  fine  grid  is  aliased 
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smooth  on  a  coarser  grid,  which  implies  that  we  do  not  want  to  move  to  the  coarse 
grid  if  any  oscillatory  components  are  present  in  the  error. 

Let  the  residual  be  defined  as 

f=f-Av,  (6.3) 

and  derive  the  residual  equation 

Ae  =  f 

by  subtracting  (6.3)  from  (6.1).  It  should  be  noted  that  relaxing  on  the  original 
equation  Au  =  f  with  an  arbitrary  initial  guess  v  is  equivalent  to  relaxing  on  the 
residual  equation  Ae  =  r  with  an  initial  guess  of  the  zero  vector. 

We  can  combine  the  above  ideas  into  a  multilevel  method  by  relaxing 
on  the  fine  grid  until  only  smooth  error  components  remain,  then  solving  the  residual 
equation  for  T on  the  coars<‘  grid,  and  finally  correcting  the  fine  grid  approximation  v 
by  that  amount.  There  are  two  advantages  of  approaching  a  problem  in  this  fashion. 
First,  we  can  address  the  troublesome  smooth  error  components  on  the  coarse  grid, 
and  more  importantly,  by  Ccisting  the  problem  on  a  coarser  grid  we  reduce  its  size  - 
in  this  example  by  75%.  In  general,  the  reduction  factor  is  ^  in  moving  from  to 
where  D  is  the  number  of  dimensions  in  the  problem.  This  is  the  general  basis 
of  coarse  grid  correction,  which  can  be  expressed  as  the  following  procedure; 

•  Relax  on  Au  =  f  for  v  on 

•  Compute  r  =  f  —  Av 

•  Solve  Ae  =  f  on 
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•  Correct  v  -f-  u  +  e  on 

Some  mechanisms  must  be  developed  for  transfering  vectors  between 
the  coarse  and  fine  grids,  as  well  as  for  representing  the  residual  equation  on  the 
coarse  grid.  This  will  require  some  additional  notation.  Superscripts  indicate  the 
grid  on  which  an  expression  is  defined,  e.g.  is  the  right-hand-side  represented 
on  the  coarse  grid,  while  is  the  approximate  solution  represented  on  the  fine  grid. 
In  addition,  the  intergrid  transfer  operators  and 

are  defined,  which  serve  to  transfer  quantities  between  the  grids.  The  operation  of 
transfering  information  from  the  coarse  to  the  fine  grid  is  called  interpolation  ,  and 
from  the  fine  to  the  coarse  grid  is  known  as  restriction  .  For  the  coarse  grid  correction 
procedure,  we  need  to  restrict  the  fine  grid  residual  to  the  coarse  grid,  interpolate  the 
coarse  grid  error  correction  to  the  fine  grid,  and  represent  the  residual  equation  on 
both  grids. 

The  interpolation  operator  is  denoted  by 

ih  ^2h  _ 
hh^  —  V  , 

and  produces  a  fine  grid  vector  whose  entries  are 

and  +  j  =  0;y-l. 

The  even  numbered  grid  points  on  are  exactly  the  grid  points,  while  the  odd 
numbered  grid  points  on  are  computed  by  averaging  adjacent  grid  points  from 
0^^,  that  is,  by  linear  interpolation.  Figure  39  shows  the  action  of 
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Interpolation  acts  on  coarse  grid  vectors  to  produce  fine  grid  vectors.  It 
is  critical  to  consider  the  smoothness  of  the  vector  being  interpolated,  as  only  smooth 
vectors  on  the  coarse  grid  can  be  accurat<‘ly  represent ^m!  on  the  fiiu'  grid.  Figure  40 
shows  the  interpolation  of  a  smooth  and  osrillator\  «Tror  vector  to  the  fine  grid.  If 
the  actual  error  is  oscillatory  on  the  fin«*  grid,  tlien  interpolation  will  not  accurately 
represent  it. 


The  restriction  operator  takos  xoctors  from  the  fine  grid  to  the  coarse 


grid  as 


There  are  several  choices  for  restriction  operators,  of  which  we  will  use  full  weight- 
which  is  defined  as 


-  4(^2i-l  +  +  ^2j+l)) 
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Figure  40.  (a)  If  the  error  (indicated  by  o  and  •)  is  smooth,  an  interpolant  of  the 
coarse  grid  error  (indicated  by  o)  should  give  a  good  representation,  (b)  If  the  error 
is  oscillatory,  an  interpolant  may  give  a  poor  representation. 


Values  for  coarse  grid  vectors  are  weighted  averages  of  values  at  neighboring  fine  grid 
points.  Figure  41  shows  the  action  of  a  very  good  example  of  why  an  oscillatory 
vector  should  not  be  restricted. 


Figure  41.  Restriction  of  a  fine  grid  vector  to  the  coarse  grid. 


Incorporating  this  new  notation,  we  formally  define  the  coarse  grid  cor¬ 
rection  scheme  as 
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•  Relax  u  times  on  A^u’^  =  on  with  initial  guess  n*. 

•  Compute  f —  A’^v^). 

•  Solve  =  r  on 

•  Correct  the  fine  grid  approximation  A 

Here,  superscripts  denote  the  grid  where  the  quantity  is  defined,  and 
the  parameter  u  represents  the  number  of  relaxation  sweeps  performed  before  moving 
to  the  coarse  grid. 

The  quantity  A^^,  the  coarse  grid  version  of  remains  to  be  defined. 
To  this  end,  assume  that  the  fine  grid  error  e*  lies  entirely  in  the  Range^I^h)-  This 
means  that  for  some  Then  we  can  write  the  residual 

e<|ualion  as 

A^c^  =  =  r'*. 

Appl  ying  the  restriction  operator  to  both  sides  yields 

i2h  ih  —2h 

‘h  •'*  ‘2i  "  =  hi  '  • 


or.  denoting  by 

which  has  the  same  form  as  the  residual  equation,  except  that  it  is  on  the  coarse  grid. 
Thus  it  is  reasonable  to  define 


=  /JM*/,**,  (6.4) 

which  is  known  as  the  Galerkin  condition  . 
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The  intergrid  transfer  operators  can  be  defined  as  matrices.  Specifically 


iltul 


-  - 
hh  -  2 


1  1 
2 
1 


1 

2 

1 


1  2  1 


1  2  1 


1  2  1 


The  intergrid  transfer  operators  satisfy  the  relationship 


cen. 


h  \T 


(6.5) 


The  relationships  (6.4)  and  (6.5)  together  are  known  as  the  variational  properties  . 
Satisfying  the  variational  properties  will  further  facilitate  the  analysis  of  the  coarse 
grid  correction  scheme. 
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Next,  consider  the  action  of  the  intergrid  transfer  operators  on  the 
modes  of  as  given  by  (6.2).  Full  weighting  applied  to  modes  results  in 

4  Zk  =  COS  {^)Zk  ,  1  <  <  y , 

and 

Here  modes  above  y  are  oscillatory,  while  those  below  ^  are  smooth.  Thus,  full 
weighting  acting  on  the  mode  of  produces  a  multiple  of  the  mode  of  A^^. 
However,  acting  on  the  {N  —  ky^  mode  of  A^  full  weighting  also  produces  a  multiple 
of  the  k‘^  mode  of  A^^  [Ref.  28].  So  full  weighting  acting  on  an  oscillatory  mode  will 
return  a  smooth  mode.  For  this  reason  it  is  essential  that  we  relax  until  only  smooth 
modes  remain  before  moving  to  the  coarse  grid. 

The  modes  of  A^^  are  given  by 


^2h 

Zk 


\ 


( sio(i«) ; 

It  can  be  shown  /citebrig87  that  the  interpolation  operator  applied  to  these  modes 
results  in 


-  sin4^)^_fc,  1  <  A;  <  ^ 


kn 


kiT 


N 


'2N 


2N‘ 
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Thus,  interpolation  produces  both  smooth  and  oscillatory  modes  on  the  fine  grid. 
Applied  to  the  smoothest  modes  on  e.g.,  modes  with  the  magnitude  of 

the  multiplier  of  the  smooth  modes  is  of  0(1),  while  that  of  the  oscillatory  modes  is 
of  0  (;^)-  Therfore,  the  result  of  interpolation  is  primarily  a  smooth  mode  on 

With  the  knowledge  of  the  action  of  the  intergrid  transfer  operators,  we 
return  to  the  coase  grid  correction  scheme,  whose  steps  are 

•  Relax  u  times  on  Cl^  using  iterative  method  P  :  v'^ 

•  Full  weight  r  ^  to 

•  Solve  the  residual  equation  exactly: 

•  Correct  the  approximation  on  PP:  n ^ 

The  process  may  be  written  as  a  single  operation,  namely 

^  Jr 

The  exact  solution  is  unaffected  by  coarse  grid  correction,  so  we  have 

Subtracting  these  two  expressions  yields 

f-  [7- A'‘]P‘'e\  (6.6) 

Denote  the  action  of  coarse  grid  correction  as  CG  :  so  that  (6.6)  can  be 

rewritten  as 

f- 
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We  wish  to  determine  the  action  of  CG  on  the  modes  of  A^.  Let  Sk  = 
sin^(^)  and  =  cos^(^).  If  we  apply  CG  without  relaxation,  it  can  be  shown 
that 

CGzk"  =  SkZk^  -\r 

and 

CG^k_k  =  ckz^k^  +  l<k<j. 

Thus  when  CG  is  applied  to  either  smooth  or  oscillatory  modes,  both  smooth  and 
oscillatory  modes  are  produced.  However,  we  must  again  look  at  the  magnitudes  of 
the  resulting  modes. 

For  k  «  A  we  have 


Thus  CG  without  relaxation  acting  on  smooth  modes  produces  both  smooth  and 
oscillatory  modes  of  small  magnitude,  while,  acting  on  oscillatory  modes,  produces 
both  smooth  and  oscillatory  modes  of  0(1).  To  prevent  using  CG  on  oscillatory 
modes,  relaxation  is  performed  first  to  eliminate  them.  Then,  after  the  fine  grid 
approximation  is  corrected,  more  relaxation  can  be  performed  to  eliminate  those 
oscillatory  modes  excited  by  interpolation. 

Now,  the  final  version  of  the  coarse  grid  correction  scheme  is 
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•  Relax  1^1  times  on  on  with  initial  guess 

•  Compute  -  A^^v^). 

•  Solve  =  r  on 

•  Correct  the  fine  grid  approximation  •(— 

•  Relax  1/2  times  on  A^u^  —  on  with  initial  guess  v^. 

The  only  unanswered  question  is  how  to  solve  the  residual  equation  on 

The  answer  is  to  think  recursively.  We  keep  transfering  the  problem  to  coarser 
and  coarser  grids  until  it  is  small  enough  to  be  solved  easily  with  a  direct  method. 
The  concept  is  illustrated  in  the  following  scheme. 

Relax  Vi  times  on  A^xi'"  =  with  initial  guess 
Compute 

Relax  ui  times  on  with  initial  guess  =  0. 

Compute  p^  = 

Solve 

Correct 

Relax  U2  times  on  A'^^iP^  =  p^  with  initial  guess 
Correct  u *  + 

Relax  1/2  times  on  A^u^  =  P  with  initial  guess  v^. 

The  scheme  telescopes  down  to  the  coarsest  grid,  which  may  be  a  single 
point,  and  works  its  way  back  up  to  the  finest  grid.  Figure  42  shows  the  schedule  of 
grids  visited  during  the  execution  of  the  algorithm  for  six  levels,  which  resembles  the 
letter  V .  For  this  reason,  the  multilevel  scheme  is  often  referred  to  as  a  V-cycle. 


105 


2h 


4h 


8h 


16h 


32h 


Figure  42.  Schedule  of  grids  for  a  V-cycle, 


To  summarize,  all  frequencies  (excepting  the  smoothest)  are  eventually 
oscillatory  on  some  grid  where  they  are  eliminated  by  relaxation.  The  smoothest 
frequencies  are  eliminated  by  the  direct  solve  on  the  coarsest  grid.  When  used  these 
schemes  are  used  together  in  the  form  of  a  V-cycle,  all  error  components  are  elimi¬ 
nated.  The  intergrid  transfer  operations  that  make  up  CG  are  chosen  to  complement 
each  other  and  the  relaxation  method  being  used.  We  want  the  restriction  opera¬ 
tor  such  that  the  smooth  components  of  the  error  on  appear  oscillatory 
transfered  to  so  that  relaxation  will  be  effective  in  eliminating  them.  Similarly, 
we  want  the  interpolation  operator  to  faithfully  represent  smooth  components  of  the 
error  on  smooth  when  they  are  transfered  to  All  of  these  operations  must 
complement  each  other,  or  the  multilevel  method  may  not  be  as  effective  as  it  could 
be. 

The  coarse  grid  correction  scheme  derived  for  this  model  differential 
equation  illustrates  the  workings  of  all  the  elements  of  a  multilevel  method.  Unfor- 
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tunately,  the  image  reconstruction  problem  is  nothing  like  the  model  problem,  and  if 
we  took  the  naive  approach  of  applying  the  elementary  multilevel  approach  outlined 
above  numerous  difficulties  would  arise.  For  example,  what  are  the  fine  and  coarse 
grids  in  this  setting?  What  do  restriction  and  interpolation  represent  in  terms  of 
the  physics  of  the  problem?  Does  restriction  make  a  smooth  error  component  appear 
more  oscillatory  in  this  setting?  Can  the  error  be  accurately  represented  on  the  coarse 
grid?  These  and  other  questions  must  be  addressed.  Clearly,  a  different  approach  is 
needed. 

To  cast  the  image  reconstruction  problem  in  a  multilevel  setting,  we 
will  use  the  abstract  multilevel  projection  method, or  PML, approach.  In  PML,  the 
problem  is  discretized  by  projections, a,nd  these  projection  operators  in  turn  define  the 
intergrid  transfer  operators  and  the  appropriate  relaxation  method. 

2.  Multilevel  Projection  Methods  (PML) 

The  multilevel  projection  methodology,  due  to  McCormick  [Ref.  2],  was 
developed  so  that  a  variety  of  problem  types,  not  limited  to  elliptic  partial  differential 
equations,  could  be  cast  in  a  multilevel  setting.  PML  is  useful  in  that  it  provides  a 
formalism  that  greatly  eases  the  development  of  intergrid  transfer  operators  and  re¬ 
laxation  schemes  that  complement  each  other.  The  designer  of  the  multilevel  scheme 
must  specify  a  set  of  subspaces,  and  the  other  components  of  the  scheme  are  deter¬ 
mined.  This  is  significant,  for  most  of  the  work  involved  in  designing  a  multilevel 
scheme  is  taken  up  in  choosing  such  components,  which  is  a  difficult  process  and 


107 


sometimes  involves  trial  and  error.  If  these  operators  are  not  matched,  the  multilevel 
scheme  might  not  be  effective. 

In  PML  the  problem  is  discretized  by  orthogonal  projections,  and  the 
projection  operators  themselves  lead  to  the  correct  choices  for  intergrid  transfer  op¬ 
erators  and  relaxation  schemes.  We  now  briefly  describe  the  general  principles  of 
PML. 

Let  H\  and  H2  be  Hilbert  spaces,  let  L  ://]—>  //2,  and  let  f  G  H2  he 
given,  and  consider  a  problem  defined  by 

Lu  =  /,  u  €  II  \. 

Define  A  (u)  =  Lu  —  f  =  0,  where  A'  :  //|  — >  Z/^.  I'heti  we  write  the  problem  as 
find  u  ^  Hi  such  that 

A'di)  =  0 

The  problem  is  generally  posed  in  tlii"  fa>hion  m»  that  it  may  be  treated  in  equation 
(strong),  variational,  or  weak  form,  although  wr  will  only  consider  equation  form 
here. 

One  of  the  basic  prinripl«*>  of  the  multilevel  projection  methodology 
is  that  discretization  is  accomplished  by  projections,  a  procedure  that  relates  the 
continuum  problem  to  a  discrete  problem  on  level  h  . 

Let  be  a  finite-dimensional  subspace  of  Z/j,  and  let  :  Hi 
be  an  orthogonal  projection  of  Hi  onto  S^.  Similarly,  let  P^''  :  H2  be  an 

orthogonal  projection  of  H2  onto  a  finite-dimensional  subspace  T^. 


108 


Assume  also  that  mappings  Psh  and  Pxh  exist  from  the  level  h  spaces 
to  the  continuum  spaces 

Pgn  :  ^  Hi,  5"*  C  Hi, 

and 

Prh  :  T*  -5-  Hi,  C  Hi. 

These  are  generally  identity  operators,  but  we  carry  them  for  generality. 

Here,  and  are  subspaces  corresponding  to  some  discretization 
parameter  h.  That  is,  they  could  be  the  space  spanned  by  a  set  of  finite  element  basis 
functions  on  a  grid  with  nodal  spacing  li:  they  could  he  continnum  functions  sampled 
on  a  grid  with  spacing  h,  or  (for  the  image  reconstruction  problem)  they  could  be 
spaces  spanned  by  strips  of  width  h. 

Finally,  assume  that  a  similar  mU  of  orthogonal  projections  exist  for 
level  2h  subspaces  and  Then  the  interlevel  transfers  betwe<*n  spaces  at  levels 
/;  and  2h  are  defined  implicitly  by  fmdiim  o|»erators  and  that  make  the  diagram 
in  Figure  43  commute. 


That  is,  and  —*  >'**  are  defined  implicitly  by 


and 


_  j2h  pS'' 


An  analogous  diagram  exists  for  the  Hi  subspaces  T*  and 
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H, 


Figure  43.  Discretization/ Coarsening  Diagram  for  Hi. 

The  process  described  above  is  known  as  discretization  by  projections^ 
and  it  dictates  what  the  coarse  grid  and  intergrid  transfer  operators  should  be.  Fol¬ 
lowing  this  approach,  the  discretized  problem  becomes 


=  0,  li  €  Pi, 


which  we  write  <is 


K^(u^)  =  0,  5^ 


It  is  important  to  note  that  is  defined  by  the  continuum  operator  I\  and  the  action 
of  the  projection  operators,  which  is  fundamental  to  PML.  That  is,  I\  . 

For  ease  of  development,  it  is  usually  assumed  that  the  subspaces  are 
conforming  ,  with  C  and  C  T^,  and  that  the  so-called  variational  proper¬ 
ties  are  satisfied,  that  is 

c€7^, 
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and 


_  j2h  L^h  jh 

^2h' 

It  is  important  that  this  discretized  problem  be  realizable, or  repre¬ 
sentable  in  a  computational  setting.  The  usual  approach  is  to  specify  the  subspaces 
5^  and  in  terms  of  finite  elements,  choose  a  basis  for  each,  and  then  rewrite  the 
discrete  problem  in  terms  of  the  coefficients  of  the  unknown  expanded  in  the  basis 
for  S'^. 

The  two  main  components  of  a  multilevel  method  needed  to  solve  this 
problem  are  relaxation  and  coarse  grid  corection.  Relaxation  will  take  the  form  of 
a  generalized  block  Gauss-Seidel  method.  To  develop  the  relaxation  method,  define 
tin*  block  subspaces  5^,  1  <  f  <  m,  such  that 

m 

/=i 

Thi.s  is  not  necessarily  a  direct  sum,  although  it  may  be.  Therefore  any  element  of 
ran  be  written  as  a  {not  necessarily  uni(iue)  linear  combination  of  the  elements  of 
S';,  e.g. 

m 

where  €  S^. 

/=i 

We  can  define  similar  block  subspaces  for  T^.  Relaxation  applied  to  the  discretized 
problem  is  then  given  by  the  steps 
For  £  =  1,2,  ...,  m 

Solve  -b  =  0. 


Ill 


Set  <r- 


Essentially,  at  each  step  we  are  seeking  an  element  from  the  appropriate  block  sub¬ 
space  such  that  after  adding  this  element  to  the  current  approximation  the  projection 
of  the  residual  onto  the  subspace  vanishes.  If  the  block  subspaces  are  chosen  to 
be  the  standard  basis  vectors,  the  relaxation  method  as  defined  here  is  just  point 
Gauss-Seidel.  Relaxation  will  be  represented  as 

The  coarse  grid  correction  procedure  is  not  difficult  to  define,  as  it 
involves  an  exact  solution  on  level  2h  .  The  procedure  is  represented  as 


ami  is  given  by  the  steps 

Solve  =  0,  ^  52/. 

Set  =  u^  + 

Here  we  seek  an  element  from  the  5^^  subspace  that  solves  the  residual  equation  on 
level  2h  .  This  element  is  then  used  to  correct  the  level  h  approximation.  We  can 
combine  coarse  grid  correction  with  relaxation  to  produce  a  two-level  PML  method 
that  will  be  denoted  as 

and  is  given  by  the  steps 
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As  before,  the  exact  solve  on  level  2h  can  be  approximated  by  a  recursive  call  4— 
PM leading  to  a  PML  V-cycle. 

B.  PML  FORMULATION  OF  THE  IMAGE  RECON¬ 
STRUCTION  PROBLEM 

We  can  apply  PML  methodology  to  the  image  reconstruction  problem  to  de¬ 
velop  a  multilevel  method  that  is  more  effective  than  relaxation  alone. 

1.  Discretization  by  Projections 

We  begin  by  showing  that  the  natural  pixel  discretization  of  Chapter 
1\’  is  in  fact  a  discretization  by  orthogonal  projections.  The  continuum  problem  to 
be  discretized  is  Au  =  /,  where  A  :  L2{^)  — >  .  We  take  //[  =  Z/2(n)  and 

II2  =  .  The  subspace  =  Range{A'‘),  the  span  of  the  characteristic  stri[) 

functions.  We  take  =  TZ^ ,  which  implies  that  =  Is-,  the  N  x  ,V  identity 

matrix. 

The  orthogonal  projection  operator  could  be  explicitly  calculated 
by  applying  the  Gram-Schmidt  process  to  the  characteristic  strip  functions,  so  that 

i=i 

where  the  ^j’s  are  a  set  of  orthogonal  basis  functions  produced  by  Gram-Schmidt  on 
the  ‘tpj's.  However,  the  following  result  shows  that  it  is  not  necessary  to  explicitly 
produce  the  orthogonal  projection  operator  P'^*. 


113 


Theorem  6.1:  For  each  angle  (j)j^  1  ^  J  ^  M ^  let  f]  be  exactly  partitioned 

into  A^i(j)  parallel f  non-overlapping  strips ^  and  let  N  =  Number 

the  strips  from  1  to  N  and  let  'il)^{x^y)  be  the  characteristic  function  of  the 
strip.  Let  be  the  subspace  of  the  Hilbert  space  Hi  =  L2{^)  spanned  by  the 
set 

Then  the  matrix  equation 

Ba  =  f 

is  a  discretization  by  orthogonal  projection  of  the  problem  Au  =  /. 

Proof:  Using  the  various  subspaces  as  defined  above,  the  discrete  equation 
will  be 

InAP^\  =  /, 

where  is  an  orthogonal  projection  of  u{x,y)  onto  5^,  and  therefore 

P^\  =  J2Qji;j{x,y)  =  A*q 

for  some  a  G  .  Since  is  an  orthogonal  projection,  we  must  have 
(u  —  P^'^u)  JL  V’j  for  every  tl'j  G  5*.  Hence  for  1  <  i  < 

0  =  (u  —  P^  u,  V’j) 

=  {u-A'a^xpj) 

N 

u  -  XI  ^k'4’kA'j 

k=\ 

N 

fc=i 

Therefore  if  P^''u  =  A'a  is  an  orthogonal  projection  of  u{x,y)  into  5^,  then 
the  vector  S  must  satisfy 


^  <  u,  01  >  '^ 

<  U,  02  > 

<  ti,03  > 


=  Bq. 


(6.7) 


<  u,0Ar  > 


UA 


But  the  left  side  of  (6.7)  is  just  equal  to  Au^  hence  we  conclude  that  a  must 
solve  Au  =  Ba^  and  the  projection-  discretized  form  of  Au  =  /  is  just 
Ba  =  /.  I 

In  this  way  our  natural  pixel  discretization  of  the  problem  can  formally 
be  considered  a  discretization  by  projection  in  the  PML  sense,  and  we  need  not  con¬ 
cern  ourselves  with  explicitly  forming  the  projection  operator.  With  this  discretiza¬ 
tion  in  hand,  the  remaining  concepts  of  PML  can  be  applied  in  a  straightforward 
manner.  We  continue  by  defining  the  coarse  subspaces  5^^  and  in  a  fashion  that 
leads  to  a  useful  multilevel  algorithm. 

2.  Intergrid  Transfer  Operators 

Let  5*  be  the  span  of  the  N  characteristic  strip  functions  where  h 
is  some  parameter  that  indicates  the  level  of  the  discretization.  For  example,  h  may 
indicate  the  width  of  the  widest  strip  function  at  that  level.  Suppose  that  there  is 
an  even  number  of  strip  functions  for  each  of  the  M  angles,  and  that  we  number  the 
functions  from  tl’x  to  i'%  in  a  way  so  that  two  adjacent  strips  on  any  view  are  always 
numbered  consecutively.  Then  the  subspace  can  be  constructed  as 

=  span  where  V’f  =  V’u-i  + 

Each  characteristic  strip  function  in  the  coarse  subspace  is  the  union  of  two  adjacent 
fine  space  strip  functions.  We  define  the  coarse  grid  problem  by  thickening  the  x-rays. 
Using  these  coarse  subspace  strip  functions  and  following  the  procedures  of  Chapter 
IV,  we  can  define  (A^^)*  :  —¥  by 
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(/t“)-a“  =  (v-;*  0“  0|‘  ■  ■  ■ 


a 


2h 


a\ 


2h 


a: 


2h 


«Ar/2  y 


2/i 


which  in  turn  defines  :  5^^  — )•  by 


= 


<  > 


<  > 


<  •4>'nI2^U  > 


An  application  of  Theorem  6.1  using  these  level  2h  subspaces  leads  to  the  projection 
discretized  coarse  level  problem  where  is  an 

(A72)  X  [N 12)  matrix  with  entries  {b]j)  = 

Having  found  orthogonal  projections  into  the  and  subspaces,  we 
next  must  determine  intergrid  transfer  operators  Ip  and  I^h  as  implicitly  defined  in 
Figure  43.  The  derivation  is  centered  around  the  definitions  of  the  coarse  subspace 
strip  functions  ‘ipp. 


Lemma  6.1:  Let  the  coarse  subspace  strip  function  ipp  be  the  union  of  two 
adjacent  fine  subspace  strip  functions  given  by  rpp  =  P'"'  1  — 
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k  <  N/2.  Then  the  operators  and  :  5^^  — >■  are  related 

by 

^2/i  _  j2h 

where  is  an  {N/2)  x  (N)  matrix  given  by 


1 


\ 


1 


Furthermore,  the  adjoint  operators  (A^)*  :  and  {A"^^)*  :  — >  5^^ 

are  related  by 

{A^^)*  =  {A^)*[ll^Y. 

Proof:  The  component  of  the  vector  A"^^u  E  is  related  to  the  {2k—  1)^* 
and  2k^^  components  of  A^u  E  by 


(V’2fc-l  +  =  (V’2Jt-n«>  +  (^2Jt.«>  =  (11) 


{i’2k^  «)  )  ' 


Now,  partition  the  vector  A^u  into  blocks  consisting  of  pairs  of  adjacent  entries 
and  form  the  matrix  by  placing  the  block  (1  1)  in  the  {2k—  1)*‘  and  (2A;)‘^ 
positions  of  the  k^^  row  of  a  {N/2)  x  {N)  matrix  of  zeros,  for  1  <  /:  <  N/2. 
The  matrix  vector  multiplication  I^^A^u  produces  A^^u  and  proves  the  first 
part  of  the  lemma. 


The  second  part  of  the  lemma  is  arrived  at  by  observing  that 


{j^2h).^2h  ^  (^^2h  ^2h  ... 

=  (V’f  +  V’2  V’s  +  V’4  ■■■  V’^/2-I  "b  V’N/2)  ^ 


-‘2h 


1 


—  ^2  V’3  ^4  ■  ■  ■  V’Ar/2-1  V’^/2) 


{A^y 


1 

1  / 
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The  second  part  of  the  lemma  verifies  that  the  operator  gives  a 
consistent  definition  to  the  adjoint  of  the  coarse  space  operator  showing  that 

=  {ii^A^y  =  {A^y[iff . 

Thus  we  can  define  I^h  —  Finally,  we  show  that  the  intergrid  transfer 

o{)erators  as  defined  satisfy  the  variational  properties. 


Theorem  6.2:  The  discrete  operators  satisfy  the  variational  properties 

4  =  d'lT- 

and 

Q2h  ^  py^B'^i!^^. 

Proof:  The  first  property  is  satisfied  with  c  =  1,  as  was  shown  in  the  proof 
of  Lemma  6.1.  For  the  second  property,  consider 

=  i}yA\A^ri!^, 

= 


I 


3.  Relaxation 

A  relaxation  scheme  can  be  developed,  following  the  principles  of  PML, 
by  partitioning  the  discrete  spaces  5"^  and  into  block  subspaces 


771 

s"  =  E5? 


£=1 


and 

i=\ 
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One  obvious  choice  of  block  subspaces  is  to  let 


S(  =  and  T/  =  span  £  =  1  :  A^. 

Here,  the  are  the  standard  unit  basis  vectors  for  .  Using  these  choices,  we  have 
the  following 


Lemma  6.2:  The  PML  relaxation  scheme  on  Au  =  /,  using  the  above  block 
subspaces,  is  implemented  by  performing  point  Gauss-Seidel  iteration  on  the 
matrix  equation  Ba  =  f. 

Proof:  The  step  of  PML  relaxation  requires  finding  a  value  s  satisfying 

pr  -  p‘)  =  0 

where  is  the  current  approximation  of  the  solution,  whicli  is  then  updated 
by 

u'  +  SV(. 

Now  xjjf  =  and  since  P"'"' u'‘  €  S''  wr  must  have  P''*' td'  =  {A'^yS^ 

for  some  €  T^.  Hence  we  need  .>«  to  satisfy 


0  =  Pi"  (AHiA'A^y-  ^  >{A^-yty;)  -  r^) 

=  Pi"  ^  .suf)  -  fl 

Realizing  that  the  action  of  the  projection  is  nothing  more  than  forming 
an  inner  product  with  we  seek  an  .»  such  that 

The  solution  is  given  by 

where  bj  is  the  row  of  B^.  Therefore,  the  step  of  PML  relaxation  is 

a  ^  a  +  -[-[}(  -bfQ  ), 

b(( 

which  is  precisely  the  correction  of  the  step  of  Gauss-Seidel  applied  to 

Bh^h  ^  jh^  I 
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4.  Coarse  Grid  Correction 


The  final  component  required  to  complete  our  PML  formulation  is  the 
coarse  grid  correction.  For  the  problem  Au  =  /,  it  is  defined  as  finding  that  element 
y2h  ^  g2h  Yvj^jch  satisfies 

-  f^)  =  0,  (6.8) 

where  is  the  current  approximation  of  the  solution  in  the  fine  space  S^.  The 
correction  is  then  given  by 

u''  + 

We  know  that 

where  represents  P^^  uK  We  also  know  that  since  €  5^^,  there  exists  a 

vector  E  such  that  Hence  AP^^^ 

Noting  also  that  then  (6.8)  becomes 

+  p2/.-2/»  _  p^up.  ^  Q 


This  establishes  that  under  the  PML  methodology,  the  coarse  grid  correction  scheme 
is  equivalent  to  that  used  for  conventional  model  problems.  In  other  words,  we  have 
proved  the  following  result. 


Lemma  6.3:  Let  =.  p-  be  the  discretization  by  projections  of  Au  =  f 

using  the  characteristic  strip  functions  ipj  and  xfjf .  Suppose  that  {A^)*a'^  is 
the  representation  of  the  current  approximation  after  relaxation.  Then  the 
PML  coarse  grid  correction  scheme  is  given  by  the  steps 
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1)  Set  =  /f  (/*  - 

2)  Solve  the  equation 

3)  Correct  the  approximation  5^  •<— 

The  complete  PML  two-level  method  is  formed  by  combining  relaxation 
and  coarse  grid  correction. 

Two-level  PML  method:  a*  i—  PML{B^,a'^,f^) 

1)  Relax  Ui  times  on  B^a^  = 

2)  Set  -  B^a'^). 

3)  Solve  the  equation 

4)  Correct  the  approximation 

5)  Relax  1/2  times  on  B^a^  = 

The  additional  1/2  relaxation  sweeps  at  step  5  are  optional,  but  have 
been  observed  to  improve  performance  in  model  problems,  so  are  included  here  for 
generality.  As  discussed  earlier,  the  exact  solve  at  step  3  can  be  replaced  by  a  recursive 
application  of  the  entire  process,  so  that  the  only  time  an  exact  solve  is  required  is 
on  the  coarsest  subspace.  To  realize  such  a  recursion  in  the  PML  setting,  define  the 
coarser  subspaces  S'-'",  for  j  =  1,2,  •••  by  taking  the  characteristic  strip  functions 
that  define  the  new  subspace  to  be  the  pairwise  joining  of  strip  functions  in  the  current 
subspace,  just  as  was  done  to  form  S^^  from  S*.  The  coarsest  level  in  this  context  is 
one  thick  x-ray  that  completely  covers  the  image  space  at  each  angle.  The  resulting 
linear  system  would  be  diagonal  and  of  size  M  x  M,  and  could  be  solved  directly.  The 
recursive  version  of  the  method,  a  PML  V-cycle,  is  given  by  the  following  algorithm: 
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PML  V-cycle:  ^  PMLV{B^,  f^) 


1)  Relax  ui  times  on  =  /^. 

2)  If  at  the  coarsest  level,  go  to  3.  Otherwise 

a) 

b)  ^  0. 

c)  <r-  PMLV{B'^^,a^^J'^^). 

d)  a*  ■<— 

3)  Relax  U2  times  on  B^a^  =  f^. 

5.  Convergence 

We  next  look  at  the  convergence  properties  of  this  multilevel  method, 
and  present  a  formal  proof.  It  is  of  limited  use  as  a  convergence  proof,  as  it  depends 
on  some  constants  that  cannot  be  determined  a  priori.  However,  the  result  has  some 
practical  applications  related  to  the  performance  of  the  multilevel  method. 

Consider  comparing  residual  norms  before  and  after  a  V-cycle  is  per¬ 
formed.  This  entails  writing  the  algorithm  in  more  detail,  so  that  residuals  can  be 
examined  at  various  steps  within  a  cycle.  It  also  involves  placing  side  conditions  on 
the  relaxation  scheme  to  measure  its  effectiveness.  Ultimately,  we  desire  the  norm  of 
the  residual  to  be  reduced  by  the  scheme  at  each  level  in  the  V-cycle. 

Define  an  artificial  level  Let 

/|  =  =  I,  B^  =  B\  and  r2  =  r\ 

2 

Note  that  this  artificial  level  is  identical  to  the  finest  grid  in  every  respect.  Then  we 
have  the  following  algorithm  for  k  levels: 
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i- 

Algorithm  MG(j,  r  2 ,  Qg  ),  where  L  =  2^~^h,  j  =  1  :  k. 

•  Compute  Tg  =  Iir2 

2 

•  Relax  i/i  times  on  =  Tq  with  initial  guess  Og 

•  Compute  rf  =  Vq  —  where  |||rf|||  <  p^|||r'2 1|| 

•  if  j  <  fc  (not  the  coarsest  grid)  then  •<— MC(j4-l,?'f ,0) 

•  Correct  •«— 

•  Compute  r2  =  Tg  —  B^a^ 

•  Relax  U2  times  on  B^o^  = 

•  Compute  Tg  =  Tg  —  B^a^  where  Ilirg  |||  <  e^|||r2  ||| 

•  Return 

The  side  conditions  occur  in  steps  3  and  8  of  the  algorithm,  and  involve 
the  constants  and  For  convergence  of  the  method,  we  require  that  these  side 
conditions  on  the  relaxation  scheme  be  satisfied.  For  the  first  relaxation,  on  the  way 
down  into  the  V-cycle,  the  factor  p^  is  the  amount  the  norm  of  the  residual  is  reduced 
after  relaxation  on  level  L  as  compared  to  the  residual  before  relaxation  on  level 
For  the  second  (optional)  relaxation,  on  the  way  up  out  of  the  V-cycle,  the  factor 
measures  residual  reduction  before  and  after  relaxation  on  level  L.  If  this  relaxation  is 
not  performed,  then  =  1.  Note  that  these  factors  are  functions  of  L,  and  change 
as  the  algorithm  moves  from  level  to  level.  We  also  want  the  restriction  operator 
satisfy  the  condition 


111(7 
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Our  restriction  operator  is  not  invertible,  so  we  approximate  it  with  its  pseudo¬ 
inverse  as 

(/f  )■'  ~  (/f 

and  require  this  approximation  to  satisfy  the  condition  instead.  Before  we  give  the 
main  theorem,  a  lemma  is  required. 

Lemma  6.4:  If  the  variational  properties  hold  with  c  =  1,  then  the  interpo¬ 
lation  operator  I21  is  norm-preserving  in  the  energy  norm. 

Proof: 


=  ll|/2V^III 


I 


The  next  theorem  follows  the  basic  outline  given  in  [Ref.  29],  and  uses  the  notation 
of  Algorithm  MG. 


ZJ.  ...  Ij 

Theorem  6.3:  Let  r  2  be  the  initial  residual  on  level  at  some  step,  where 
Lj  =  1  <  j  <  kf  where  k  denotes  the  number  of  levels  of  the  scheme. 

Let  (/^)“*  :  be  approximated  by  its  pseudoinverse  Further^ 

assume  that  there  exists  a  5^  £  TZ  such  that 

|||(/-(/“)'/f)„l||<J'-|||u|||,  U€SI‘-. 

Define  Ei  =  e^ft^  and  Ej  =  +  Ej^i),  j  >  1,  where  L  =  2^~‘h. 

rAenllK/iMlIISB.IIIr^lll. 

2 


Proof:  The  proof  is  by  induction  on  j.  For  j  =  1,  we  have 


J)'''3‘lll  =  lll'■3‘lll  <  f'Mlr: 
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but  since  there  is  only  one  level,  =  Tj  .  So 

IIK/ijMiii  <  £‘iiKiii  <  tV^IlKlII  =  B,|i|r«iii. 

2 

Now,  assume  the  claim  holds  for  all  levels  j  <  k.  Then 

IIK^MIII  =  \wir-m  = 

L 

since  the  operator  II  is  norm  preserving.  Therefore 

lll(^i)'’-3''lll  =  jllkilll  <  5^''ll|r2''lll 

=  -  Uiy  ll'-yi  + 

<  h'-IIK/  -  Vf)' IV-)A\\\  +  III 

<  5<'-i'-|ll-'Mll  +  ie'-|||(/f)'r|'-||| 

<  5<'-'5'-|lkfl||  +  5«'-£,-,ll|rf||| 

=  5£V''  +  £3-.)llk!'lll 

<  5fV('5''  +  £,-i)llk*lll 

=  jEjIlkhll 

I 

The  constants  e^,  5^  and  Ej  indicate  the  performance  of  the  multilevel  routine. 

As  long  as  Ej  <  1,  the  method  is  converging.  Unfortunately,  there  is  no  a  priori  way 
to  determine  the  values  of  these  constants,  as  they  change  from  level  to  level  and 
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cycle  to  cycle  as  the  algorithm  is  executing.  They  can,  however,  be  calculated  and 
monitored  during  execution  so  that  some  idea  of  the  performance  of  the  method  can 
be  gained. 

C.  ANALYSIS  OF  PERFORMANCE 

We  compare  the  performance  of  Gauss-Seidel  alone  and  the  PMLV  method 
by  considering  work  required  to  reduce  the  norm  of  the  residual.  Let  one  sweep  of 
Gauss  Seidel  on  the  finest  level  be  defined  as  one  work  unit  ( WU),  which  is  an  0{N^) 
operation.  The  work  required  for  one  V-cycle  can  be  computed  in  the  following 
fashion.  At  each  level,  we  perform  U\  +  1^2  sweeps  of  Gauss-Seitlel  and  compute  a 
residual.  Since  computing  the  residual  is  an  O(.V^)  operation,  we  lei  it  be  equivalent 
to  the  work  of  one  sweep  of  Gauss-Seidel.  ;\s  the  problem  is  roars(Mied,  the  size  of 
the  matrix  equation  to  be  solved  is  red\i<  <Ml  by  a  factor  of  *1  at  each  level.  So  for  a 
\’-cycle,  the  work  required  is 

(i^i  +  1/2  +  I )(  I  ^  77  ^  •  •  *). 

I  If* 

For  these  tests  we  use  i/j  =  2  iterations  ^^oing  into  the  V-cycle  and  1/2  =  1 

iteration  coming  out  of  it,  so  one  \^-cych*  retjulres  approximately  y  WU.  In  all  cases 
the  problem  is  coarsened  to  the  coars(‘sl  possible  h'vel,  i.e.  one  ray  per  view  for  a 
problem  of  size  A/  x  M  at  the  coarsest  level.  Figure  44  compares  the  performance  of 
Gauss-Seidel  to  the  PMLV  algorithm  for  a  problem  of  geometry  32  detectors  over  20 
angles.  This  performance  is  typical  of  that  obtained  from  numerous  experiments. 
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Porformanc®  Plots 


Figure  44.  Comparison  of  the  performance  of  Gauss-Seidel  and  PMLV. 

It  is  clear  from  the  figure  that  PMLV  initially  performs  superior  to  Gauss- 
Seidel,  but  then  this  increase  in  performance  stalls  out,  so  that  eventually  both  rou¬ 
tines  perform  in  about  the  same  manner.  However,  examination  of  the  slopes  of  the 
curves  indicate  that  further  iteration  may  favor  PMLV.  We  return  to  the  singular 
value  decomposition  to  analyze  this  behavior.  Consider  the  problem  Ba  =  /,  whose 
exact  solution  in  the  leaist  squares  sense  is 

Q  =  B'/, 


where  is  the  psuedo-inverse  of  B  .  In  terms  of  the  SVD,  a  can  be  expressed  as 


a  = 


>  —ViUi 


r  ^  -  7* 

/  =  2^— Vi  <  Ui,f  >, 

i^i 


where  r  =  rank(5).  If  there  is  measurement  noise  in  the  data  we  are  given  to  recon¬ 
struct,  so  that  instead  of  /  we  have  /  +  £,  then  solution  components  corresponding  to 
small  singular  values  will  magnify  this  noise.  Problems  of  this  nature  are  referred  to 
as  ill-posed  .  Components  in  the  near  null  space,  i.e.,  those  with  small  singular  values. 
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are  those  components  that  are  slow  to  to  be  recovered.  Thus  continued  iteration  after 
the  procedure  stalls  in  an  attempt  to  recover  these  slow  components  has  the  potential 
to  corrupt  the  solution  with  magnified  noise  [Ref.  30,  31]. 

Such  problems  require  some  form  of  regularization  to  prevent  the  ill-posedness 
from  completely  corrupting  the  approximation.  One  way  to  regularize  the  problem 
[Ref.  14]  is  to  simply  stop  iterating  when  the  algorithm  begins  to  stall.  An  ad 
hoc  approach  to  this  is  to  measure  the  difference  between  successive  residual  norms, 
and  stop  iterating  when  a  tolerance  is  achieved.  Perhaps  a  better  stopping  criterion 
exists.  Recall  from  Theorem  6.3  that  the  PMLV  algorithm  reduces  the  norm  of 
the  residual  at  each  step  by  some  factor  Ej.  This  number  can  be  computed  as  the 
algorithm  is  executing,  and  can  then  be  used  to  monitor  its  performance.  It  has 
b«*«*n  experimentally  observed  that  when  the  algorithm  begins  to  stall,  Ej  becomes 
larger  than  one  in  magnitude.  It  could  be  postulated  that  a  stopping  criterion  for 
the  iteration  is  to  monitor  the  magnitude  of  Ej,  and  then  terminate  execution  when 
Ej  >  1.  For  the  above  example,  Ej  became  greater  than  one  in  magnitude  after  the 
third  \'-cycIe,  which  coincides  with  the  stalling  of  convergence  in  Figure  44. 

The  PMLV  algorithm  applied  to  the  natural  pixel  discretized  problem  recon¬ 
structs  images  quite  well.  The  following  series  of  figures  depicts  actual  and  recon¬ 
structed  images  for  two  brain  phantoms  and  a  woman’s  face.  Note  that  this  last  image 
is  not  the  type  of  image  for  which  the  method  is  developed,  but  is  included  so  that 
the  reader  unfamiliar  with  radiological  imagery  has  something  familiar  to  observe. 
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PMLV  Conv«rg«nc«  Constant 


Figure  45.  Plot  of  the  PMLV  convergence  constant  Ej  by  sweep. 

In  the  next  chapter,  a  multilevel  fast  adaptive  composite  method  (FAC)  ap¬ 
proach  will  be  investigated  and  applied  to  the  spotlight  CAT  problem,  which  involves 
getting  high  resolution  in  one  piece  of  a  larger  image  without  discretizing  the  global 
[iroblcm  to  that  level  of  resolution.  The  natural  pixel  discretization  approach  will  be 
taken  throughout  this  development  as  well. 
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Figure  47.  Actual  and  Reconstructed  Images  -  Brain  Phantom  2. 


Figure  48.  Actual  and  Reconstructed  Images  -  Face. 
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VII. 


SPOTLIGHT  COMPUTED 
TOMOGRAPHY 


The  image  reconstruction  techniques  investigated  to  this  point  reconstruct  the 
entire  region  of  interest  over  which  the  data  has  been  collected.  In  many  cases,  more 
detailed  information  may  be  required  over  a  particular  sub-area  of  the  region.  For 
example,  a  tumor  might  be  suspected  and  the  doctor  wants  a  closer  look.  Collec¬ 
tion  of  data  over  only  the  region  containing  the  possible  tumor  result  in  inaccurate 
images,  however,  as  only  the  sub-area  is  fully  scanned.  To  overcome  this  problem, 
a  multiresolution  technique  known  as  spotlight  computed  tomography  (CT)  can  be 
employed.  [Ref.  32,  33] 

To  utilize  the  spotlight  technique,  the  sub-area  is  x-rayed  at  high  resolution, 
while  the  remainder  of  the  image  is  x-rayed  at  a  lower  resolution.  The  collection 
of  high  resolution  data  over  only  the  sub-area  reduces  the  size  of  the  resulting  linear 
system  dramatically  compared  to  uniform  high  resolution  discretization.  This  reduced 
size  in  turn  allows  the  problem  to  be  solved  in  less  time  with  fewer  resources. 

A.  NATURAL  PIXEL  DISCRETIZATION 

We  again  take  the  natural  pixel  approach  to  discretizing  the  problem.  For 
ease  of  development,  we  initially  restrict  the  refinement  to  two  levels,  and  consider  the 
general  case  later.  Let  u{x,  y)  be  the  density  function  of  the  image  to  be  reconstructed. 
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defined  in  a  square  region  Q,  of  unit  area.  Let  the  sub-area  to  be  spotlighted  be 
defined  in  a  square  refinement  region,  Q,r,  that  is  contained  within  Cl.  Assume  that, 
following  the  methodology  of  Chapter  IV,  the  image  is  x-rayed  at  coarse  resolution 
over  M  angles  with  Ni{j)  source/detector  pairs  per  angle,  1  <  J  <  M,  for  a  total  of 
^  ^lO)  x-rays  that  completely  cover  the  image  at  each  angle. 

Refinement  is  achieved  by  x-raying  Q,r  at  a  fine  resolution.  (However,  it  is 
assumed  that  all  the  data  is  collected  at  once.  The  x-ray  data  for  the  fine  grid  is 
not  acquired  at  a  later  time).  We  also  assume  that  these  fine  rays  exactly  partition 
the  coarse  rays  over  Qr,  and  that  there  are  ^2(7)  such  rays  per  angle,  for  a  total 
of  P  =  A^2(i)  fin‘*  rays.  Therefore  the  fine  rays  completely  cover  Qr  at  each 

angle.  The  following  sijnple  example  will  serve  to  illustrate  this  concept.  Assume  a 
geometry  of  three  coars<'  rays  over  each  of  two  angles,  and  that  the  refinement  region 
is  contained  within  the  inters<Ttion  of  the  center  rays  of  each  view,  as  shown  in  Figure 
49  (a).  Refinement  is  accomplished  by  dividing  the  center  rays  in  half,  as  shown  in 
Figure  49  (b). 

Note  that  the  coarse  rays  are  global  in  nature,  completely  covering  the  image,  while 
the  fine  rays  are  local,  completely  covering  only  Qr. 

Once  again,  let  the  ray  paths  be  thought  of  as  natural  pixels,  and  introduce 
characteristic  strip  functions  corresponding  to  these  pixels.  Let  for  j  = 

1  ;  A^,  be  the  coarse  strip  function,  and  i/’^,  k  =  1  ;  P,  be  the  the  k*^  fine 
(refinement)  strip  function.  Here  we  let  the  superscripts  indicate  which  resolution  is 
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(a)  (b) 

Figure  49.  An  example  of  grid  refinement. 


being  considered,  i.e.  h  denotes  fine  strips,  2h  denotes  coarse  strips,  and  let  h  denote 
a  composite  combination  of  both  resolutions  together. 

Define  the  operator  A-  :  H  7^^+^  by 


<  > 

<  u  > 


A^ 


<  ^  > 

<  > 

<  > 


A2h^2h 

A^u'^ 


^  <  V’p)  ^  >  j 


where  H  7^^  and  A^  :  H  ^  11^. 
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The  data  vectors  are  modeled  by  the  integral  operators 


u{x,y)il)]^dxdy  =  \  <  j  <  N, 

and 

J^^u{x,y)ip^dxdy  =  fk,  I  <  k  <  P, 

giving  the  system 


u(x,  y)xl)j^dxdy 

Jtz^ 

u{x,y)'il;f'dxdy 

r2h 

J2 

= 

In^ 

u(x,  y)'il)^dxdy 

!tv 

u{x^  y)'il)^dxdy 

n 

u{x^  y)xp2dxdy 

fi 

{  In'^ 

u{x^y)rppdxdy  j 

^  fp  J 

where  the  and  fj^  are  the  coarse  and  fine  ray  projection  data  for  the  image, 
respectively.  Following  the  course  of  action  outlined  in  Chapter  IV,  we  seek  the  least 
squares  solution  to  this  system,  which  is  given  by  (/l-)*a-,  where  a-  €  defined 

as 


Q—  = 


.  \ 

S2h 


a'^  j 


134 


solves  the  system  A-{A-ya-  =  Once  a-  is  found,  the  least  squares  solution  is 
given  by 

u{x,y)  = 

The  operator  [A-)*  :  —^His  defined  by 


2.^ 


{A^ya^  =  •••  V’N  i’l  ^>2  •••  V’p] 


ar 


2h 


a 


^2h 

Qn 


ai 


\^P  / 


-.A 


a 


V  / 

=  EofV-f + 

J=1  k=\ 


where  {A^^}'  :  TV^  — >  H  and  {A^y  :  Rf  — >•  H  are  defined  in  the  expected  way. 

Thus  the  image  density  is  represented  as  a  linear  combination  of  the  charac¬ 
teristic  strip  functions  by 


«(^>2/)  =  I] 

i=l  A:=l 


(7.1) 


135 


Define 


Bk  = 

where  B-  :  .  The  entries  in  B-  can  be  calculated  by  substituting 

(7.1)  for  u{x,y)  in  A-u-  =  yielding 

+  =  /f’  i  =  1:A' 

i=l  £=\ 

/,(l;a?V?‘  +  l;Q?<AfW  =  ft,  <■■  =  l:P, 

i.l  fel 

which  can  be  expanded  as 


i=l  ^=1 

finally  yielding 


=  //'••  ./  =  1  : 


=  /:•.  A-  =  1  :  r. 


1=1  ^=1 

>  =  /^  ^-  =  1  : 

1=1  ^=1 

This  last  expression  is  a  block  linear  sysl<*in,  whirh  ran  be  written  as 


ryh  -*/i 

t>-a-  = 


Q2h2h  Q2Kh 
Qh2h 


^,k\ 


/ 


O 


) 


f2H  ^ 
/ 


Thus  B-a-  =  /-  is  the  natural  pixel  discretization  of  the  spotlight  CT 

problem  A-u-  =  /^,  where  u-  =  {A-)*aK 
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It  is  possible  to  refine  the  fine  strip  functions  again  in  the  same  fashion,  which 
allows  for  recursive  refinement  to  as  fine  a  level  of  discretization  as  is  needed  to  resolve 
the  image.  Before  investigating  how  to  best  solve  this  linear  system,  we  first  analyze 
the  matrix  B-,  characterize  its  null  space,  and  discuss  some  of  its  other  interesting 
properties. 


B.  PROPERTIES  OF  THE  SYSTEM  MATRIX 

In  our  analysis  of  the  matrix  much  of  the  theory  previously  developed  in 
Chapter  IV  will  be  directly  applicable.  The  following  is  one  such  result,  the  proof  of 
which  follows  directly  those  of  Lemma  -1.1  and  Tht'orem  1.2. 

Lemma  7.1:  The  matrix  B^  is  non-tKffa/irt .  >yrnttn  Inc.  aud  positive  semi- 
dcfinite. 


From  above,  we  know  that  B-  can  be  written  io. 


B^  = 


\ 

HIKIK  IfiKK 


ip.u.  IfU. 


and  that  each  of  these  four  blocks  has  a  bliM  k  strtirture  of  its  own.  The  matrix  B^^^^ 
is  formed  from  the  intersections  of  the  coarse  rays  with  themselves,  and  is  exactly 
the  matrix  B  analyzed  in  Chapter  IV.  B'‘'‘  is  formed  from  the  intersection  of  the  fine 
rays  with  themselves.  Like  it  is  block  M  x  M  with  the  diagonal  blocks  being 

diagonal  matrices.  However,  it  lacks  the  summability  properties  of  pjnally, 

the  off-diagonal  block  matrices  B^^^  =  are  formed  by  the  intersections  of 
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Soodiqm  CT  .Vlainx 


^xxn\nxx%i 

■m^vwxx^i 


rtz  •  30872 


! 


Figure  50.  Block  structure  of  the  Spotlight  matrix. 
coarse  rays  with  fine  rays.  For  example,  would  have  the  form 


B 


hh  _ 


d/i/i 

^11 


ryfib, 

D21 


jDhh 


TDhk 

^12 


p/i/i 

^22 


ryhk 

^M2 


p/i/i 

^IM 


TDhh 

^2M 


MM 


where  B^J^  has  elements  corresponding  to  the  areas  of  intersection  of  the  fine  rays  at 
angle  4>i  with  the  fine  rays  at  angle  4>j.  The  other  two  blocks  have  a  similar  structure 
and  similar  interpretation.  The  sizes  of  the  blocks  are  is  Ni{i)  x  Ni{j),  B^h  is 

N2{i)  X  A^2(i))  is  Ni{i)  x  A^2(i))  a^nd  Bf^^  is  N2{i)  x  Ni{j).  Figure  50  illustrates 
the  block  structure  of  a  typical  matrix,  in  which  only  the  nonzero  entries  appear  in 
black. 
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The  off-diagonal  block  matrices  possess  summability  as  given  by 


Lemma  7.2;  The  off-diagonal  blocks  of  exhibit  the  following  summability 
properties: 

a)  Let  Tjt  be  the  row  of  Bh.  The  elements  of  r^  in  any  block  Bf^^  of  B^^^ 
sum  to  the  value  of  the  corresponding  diagonal  element  in  the  row  of 

Qhh 

b)  Let  fjt  be  the  k*^  row  of  B^.  The  sum  of  the  elements  of  rk  in  any  block 

Bfh^  of  B^^^  is  equal  to  the  sum  of  the  elements  in  r^  in  block  of 

Q2h2h  corresponding  to  the  coarse  rays  from  which  they  were  partitioned. 

c)  Let  ffc  be  the  k*^  row  of  5-.  The  sum  of  the  elements  of  rk  in  any  block 
B^j'^  of  B^^^  corresponding  to  coarse  rays  that  are  refined  is  equal  to  the 
sum  of  the  elements  in  rk  in  block  Bfh  of  B'^^ . 

d)  The  above  three  results  hold  for  the  columns  as  well. 

Proof:  To  prove  part  a,  consider  how  the  elements  of  row  k  for  block  Bjj^ 
corresponding  to  angle  are  formed.  The  entries  of  row  k  for  this  block  are 
the  areas  of  intersection  of  fine  strip  k  with  all  of  the  coarse  strips  for  angle 
Oj.  Since  the  coarse  rays  at  angle  (f)j  must  cover  the  entire  image,  they  must 
cover  fine  strip  k  as  well.  To  prove  part  6,  note  that  the  elements  of  rk  in 
are  the  areas  of  intersection  of  the  fine  strips  at  angle  (f>j  with  coarse  strip  k 
from  angle  (f),.  The  elements  of  r*.  in  52/12/1  corresponding  to  the  coarse  rays 
at  angle  <j)j  which  are  refined  are  their  areas  of  intersection  with  coarse  ray  k 
at  angle  </»,.  Since  the  fine  rays  are  an  exact  partition  of  these  coarse  rays,  the 
sums  must  be  equal.  This  is  geometrically  illustrated  in  Figure  51.  For  part 
c.  again  note  that  the  elements  of  r*  in  B^^^'  corresponding  to  coarse  rays  at 
angle  that  are  refined  are  their  aretis  of  intersection  with  fine  ray  k  from 
angle  4>i‘  The  elements  of  rk  in  5,^^  are  the  areas  of  intersection  of  the  fine 
strips  at  angle  <l>j  with  fine  strip  k  at  angle  <p,.  Due  to  exact  partitioning,  these 
sums  must  be  equal.  This  is  geometrically  illustrated  in  Figure  52.  Part  d 
follows  from  symmetry.  | 


Another  summability  property  exists,  based  on  how  the  refined  rays  are  de¬ 
fined,  that  has  a  role  in  determining  the  rank  of  B-.  Recall  that  by  assumption  the 
fine  rays  exactly  partition  the  coarse  rays.  Refinement  is  carried  out  by  selecting  a 
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2h2h  2hh 

B 

Figure  51.  Geometry  of  Lemma  7.2(b). 


phh 

Figure  52.  Geometry  of  Lemma  7.2(c). 

coarse  ray  and  dividing  it  into  a  number  of  thinner  rays.  Now,  consider  the  row 
of  the  matrix  B-.  The  entries  in  this  row  are  determined  by  finding  the  areas  of 
intersection  of  the  ray  path  with  all  +  P  ray  paths  defined  in  the  geometry  of 
the  problem,  taken  in  sequence.  This  leads  to  the  following  theorem. 


Theorem  7.1:  The  sum  of  the  rows  of  B-  corresponding  to  a  set  of  fine 
rays  that  were  formed  by  subdividing  any  one  coarse  ray  equals  that  row  of  B- 
corresponding  to  the  coarse  ray  in  question. 
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Proof:  Let  be  a  coarse  strip  function  which  is  refined  into  r  fine  strip 
functions,  such  that 

i+r 

V’f  =  E  V’f. 

»=j+i 

Let  be  the  row  of  B-,  given  by 


J*h 

h 


T 

<E;i;+,  </’?,<>  \ 

<  > 

<  E!;;+,  > 

<  ,  V’N  > 

<  Eiiit,  </’?>!?  > 

<  > 

<  es;+,  > 

<  V’f ,  V’2  > 

<  ESI+i  v-f,</’5  > 

V  <  V’f ,  V'p  >  y 

Expanding  the  sums  yields 


<  V’i+nV’l'"  >  +  <  >  + 


’,h 


%^>-\-<  4’-+2,  >  + 

^  ^  _L  ^  Jyh  ./,/l 


<  >  +  <  V’j>2>V'f 

<  ^i+nV’^  >  +  <  ^i+2’^^ 


>  + 
>  + 


<  V’j+i ,  V’p  >  +  <  V'j+2»  ^'p  >  + 

'j+2  + 


%1 + 


^+r- 


+  <  V’t+r,^? 
+  < 


2/.  >  \ 
> 


+  <  V’ 


+  <  > 
+  <  0jV’^2  > 


+  <V’E’^’P>  / 


An  immediate  implication  of  this  result  is  the  following 

Corollary  7.1:  The  rank  of  can  be  no  greater  than  A^  +  P  less  the  number 
of  coarse  rays  that  are  refined. 

Proof:  Let  m  be  the  number  of  coarse  rays  that  are  refined.  Each  of  these  m 
rays  leads  to  a  dependent  set  of  rows  in  the  matrix  P-,  as  per  Theorem  7.1. 
Therfore  the  rank  of  P-  cannot  exceed  +  P  —  m.  | 
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These  summability  properties  and  some  additional  analysis  will  allow  us  to 
determine  the  rank  of  B-  and  characterize  its  null  space.  Assume  v-  €  NS{B-). 


Then 


B^^  = 


Q2h2h  Q2hh 
^h2h  Qhh 


\ 

( \ 

V 

y 

y  V  J 

=  0. 


Matrix  multiplication  results  in 


Q2h2h;^2h  _|_  Q2hh;^h 

=  0 

(7.2) 

h2  hi  ^2h  1  h  h  ~*h 

D  V  -t  D  V 

=  0 

(7.3) 

We  now  investigate  under  what  conditions  (7.2)  and  (7.3)  are  satisfied.  We 
saw  from  Theorem  4.3  that  a  vector  is  in  the  null  space  of  the  matrix  B  = 
if  and  only  if  it  is  constant  by  angle  with  the  constants  summing  to  zero.  That  is, 
the  image  at  each  view  is  a  shade  of  grey,  and  when  all  views  are  superimposed  the 
result  is  a  black,  or  invisible,  image.  The  concept  of  a  vector  being  constant  by  angle 
liolds  for  the  composite  matrix  B-  as  well,  with  some  modifications.  Consider  the 
characteristic  strip  functions  for  a  refined  image  over  one  angle,  as  shown  in  Figure 
53. 


A  vector  v-,  where 


+ 


_ I  I _ I  h 

h  2h 

Figure  53.  Characteristic  strip  functions  for  a  refined  image  over  one  angle. 
is  constant  by  angle  with  respect  to  the  matrix  correspnding  to  this  geometry  if 

V-  =  (ofi  0^2  Qf3  7i  7i  72  72)^, 

where 

o,  =  02+71  =  03  +  72- 

The  contribution  in  the  refinement  region  from  a  coarse  strip,  together  with  the 
contributions  from  its  corresponding  fine  strips,  equals  the  contributions  from  those 
coarse  strips  not  in  the  refinement  region.  The  overall  result  is  that  the  composite 
collection  of  strips  forms  a  uniform  grey  image. 

We  can  generalize  this  idea  for  a  geometry  of  M  angles  by  requiring  the  com¬ 
posite  subimage  at  each  angle  to  be  uniform  grey.  Define  aj  to  the  constant  for  the 
unrefined  strips  from  angle  (fj,  ajk  for  k  =  1  :  ruj  to  be  the  values  of  the  mj  coarse 
strips  that  cover  the  refinement  region  from  angle  (j)j,  and  'yjk  for  k  =  1  :  mj  the 
values  of  the  refinement  strips  at  angle  (f)j  that  partition  the  k^^  coarse  strip  in  the 
refinement  region. 
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Definition  7.1:  Let  v-  =  where 

=  (iBi  zI52  •  •  •  wmV  and  v'*  =  {yi  y2  ■  ■  ■  mV ■ 

Subvector  Wi  is  of  length  Ni(i)  and  subvector  yj  is  of  length  A^2(i)-  Let  m 
be  the  total  number  of  coarse  strips  that  are  partitioned,  and  let  nij  be  the 
number  of  such  strips  at  angle  <f)j,  so  that  Then  v-  is  composite 

constant  by  angle  with  respect  to  B-  generated  over  M  angles  if 

■Wj  =  {aj  -  •  •  aj  ocji  dj2  •  •  •  otj  ocjV 

and 

Vj  —  {iji  ■■■  Til  7i2  ■■■  7i2  limj)  ? 

subject  to  the  constraints 

Otj  =  dtji  +  'Yji,  i  —  \  nij. 


Obviously,  if  v-  is  composite  constant  by  angle,  then  the  image  it  defines  is 
a  constant  image  with  a  value  equal  to  the  sum  of  the  q,’s  which  define  the  M 
uniform  grey  subimages.  Vectors  in  the  NS{B-)  are  composite  constant  by  angle 
and  correspond  to  constant  images,  as  will  be  shown  in  the  following  theorem. 


Theorem  7.2:  v-  €  NS{B-)  if  and  only  if  it  is  composite  constarit  by  angle 
with  ZiL  =  0- 

Proof:  Let  u-  be  composite  constant  by  angle  with  0.  VVe  can 

write 


B^^  = 


Q2h2h^h  ^  Q2hh;^h 
Qh2h^h  ^  ghh-h 


£2/i 


Consider  the  contribution  of  B-v-  corresponding  to  angle  (j)j  toward  the 
component  of  This  can  be  expressed  as 

ruj 

^  ^jkPit(k)  +  Y^in^ik!'  +  •  •  •  +  IZ 

k^SV^^  fc=l  k\  km. 

r(jt)€n'* 

which  after  applying  Lemma  7.2(b) 
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=  o.  E  E  E 

/ken^*  fc=i  *=i 

f(fc)en'>  i(k)£Q.>' 


2h2h 

'it(k) 


rrij 

=  E  /5““  +  E  0flH) 

jken^h  fc=i 

f(jk)Gn'» 

—  1 


by  Lemma  4.2(b)  and  the  fact  that  v-  is  composite  constant  by  angle.  Since 
angle  <j)j  was  arbitrarly  chosen,  we  have 


—  Aj(^l  +  0^2  +  •  •  •  +  O'M)  —  0- 

Since  row  i  was  arbitrarily  chosen  as  well,  we  have  =  0. 

Now,  consider  the  contribution  of  B-v-  corresponding  to  angle  4>j  toward  the 
component  of  i*.  This  can  be  expressed  as 


E  +  E  + 

,=1 


+  '^'KjmjPkl 


which  after  applying  Lemma  7.2(c) 


=  a,  E  +  E  +  E 

,gn2h  1=1  >=i 


—  E  E 

•=! 

/(i)6n'’ 


=  fikkOtj, 


by  Lemma  7.2(a)  and  the  fact  that  v-  is  composite  constant  by  angle.  Since 
angle  <f)j  was  arbitrarly  chosen,  we  have 

zl  =  ^kk{oc\+Oi2+  •••  Tom)  =  0. 

Since  row  k  was  arbitrarily  chosen  as  well,  we  have  i*  =  0.  Therefore, 

Now,  assume  v-  €  NS{B-).  Then  v-  €  NS{{A-)*)  and 

=  Et-?Vf +  E«>*  =  0. 

1=1  j=\ 
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h 

k-^2 


Figure  54.  Geometric  illustration  of  Theorem  7.2. 


Consider  three  adjacent  strips  and  from  the  same  projection 

angle.  Now,  select  points  (xi,t/i)  G  Sl^nCli  and  {x2,y2)  €  where  f2i 

is  the  intersection  of  (M  —  1)  coarse  grid  strips,  one  from  each  profile  and  none 
from  the  profile  containing  Sl^  and  and  is  not  a  part  of  the  refinement 

region.  This  selection  can  always  be  made  (see  Figure  32).  Therfore  we  have 

=  '^v^U’f^ix2,y2H'^v^rlj{x2,y2)  =  0. 

1=1  j=l  i=l  }=l 


Since  the  V’,-  and  tfj  are  characteristic  functions,  we  can  write 

+  +  =  0, 

ieii,  icfli 


(7.4) 


whicli  implies  that  Therefore,  outside  the  refinement  region  the 

coarse  strips  are  constant  across  angles. 

Now,  select  a  point  (xa,  1/3)  such  that  (xj.yz)  €  D  Q2  and  (x3,y3)  € 

C  ^2,  where  Q2  is  the  intersection  of  (Af  —  1)  coarse  grid  strips,  one 
from  each  profile  and  none  from  the  profile  containing  and  5'*+2-  This 
selection  can  always  be  made.  Superimpose  all  of  the  strips  over  the  image 
at  once,  forming  a  grid  of  polygons.  Since  each  coarse  grid  profile  completely 
covers  the  image,  a  point  in  the  interior  of  any  polygon  is  contained  in  a  strip 
from  each  of  the  M  profiles.  A  point  on  the  edge  (not  a  vertex)  of  a  polygon 
that  separates  a  coarse  grid  and  fine  grid  strip  will  be  contained  in  ^2?  with 
the  edge  separating  strips  and  5jt+2-  Moving  a  distance  e  to  either  side 
and  perpendicular  to  the  edge  will  locate  points  (x2,y2)  and  (x3,y3).  This 
geometry  is  illustrated  in  Figure  54. 
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Therefore  we  have 


^i;?V.^'^(2:2,y2)+^UyV’i(2;2,J/2)  =  '"V’.- '"(^3,  ^3)+$^  V’i  (2:3,  J/s)  =  0. 

i=l  j=l  t=l  j=l 

Since  the  and  are  characteristic  functions,  we  can  write 

^  vf^  +  ^  Vi^  +  +  ^fc+2  = 

which  implies  that  vl^  =  +  ^fc+2-  repeat  this  argument 

for  the  next  fine  strip  that  is  adjacent  to  5'fc+2)  finding  that  its  value  added 
to  the  value  of  the  coarse  strip  it  partitions  must  also  equal  Proceeding 
in  this  fashion  across  all  fine  strips  for  the  selected  angle  will  result  in  a  view 
whose  composite  strips  functions  all  have  the  same  value  of  Since  the 
angle  was  arbitrarily  selected,  v-  is  composite  constant  by  angle.  From  (7.4) 
we  have  Xligni  =  0.  Since  there  is  a  from  <’afli  of  the  M  profile 

contained  in  this  expression,  the  constants  must  sum  to  zero.  | 


All  immediate  consequence  of  this  theorem  is  that  image's  in  the  NS{B-)  are 
invisible.  Another  consequence  is  the  folhming  tli*i>rem.  which  relates  the  rank  of  B- 
to  the  geometry  used  to  x-ray  the  image. 

Theorem  7.3:  Let  B-  €  bf  thf  rovipositt  natural  pijd  discretized 

matrix  formed  at  M  angles,  and  a.v'uno  that  tin  njitnnnnt  rtgion  is  covered 
by  a  total  of  m  coarse  strip  funrt\on>  }‘ubdirtdtd  into  fine  strip  functions  as 
outlined  in  the  above  discussion.  Tinn  tbt  rank  of  B-  ts  .\-hB-(.\f+m-] ). 

Proof:  Consider  the  degrees  of  fret'elimi  in  s«*l«*<ting  the  values  for  the  strip 
functions.  For  the  views,  the  first  {.M  -  I )  \“alm*s  can  be  arbitrarily  chosen, 
after  which  the  final  value  is  determine«l  m»  that  the  M  values  sum  to  zero. 
The  values  of  the  m  coarse  strip  functimis  which  are  refined  are  arbitrary  as 
well,  for  a  total  of  M  -|-  m  —  1  degr«'es  of  fretHlom.  Hence  the  rank  of  B-  is 
-b  P  -  (M -1- m  -  1).  I 

The  above  results  can  be  used  to  construct  a  basis  for  NS{B-).  A  null  space 
vector,  in  terms  of  Definition  7.1,  must  satisfy 

Wj  =  (oj  •  •  •  exj  dij\  Q!j2  •  •  •  ^jrrij  OLj  •  •  •  Q!j) 
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and 


Vj  (7jl  *  *  *  7il  7j2  *  *  *  7j2  ‘  *  *  Ijmj  ’  *  ’  7jmj)  ? 
subject  to  the  constraints 


Otj  —  OZ-jx  "h  ^  —  1  • 


and 

M 

E“j  =  0- 

i=i 

Given  this  form  of  a  null  space  vector,  the  basis  is  constructed  as  follows.  The 
first  M  —  1  basis  vectors  are  obtained  by  setting  xP’^  equal  to  the  M  —  \  basis  vectors 
of  B,  while  maintaining  =  0.  This  has  the  effect  of  letting  all  the  jSjk  =  0  and 
all  the  Ojk  =  Oj-  To  construct  the  remaining  m  basis  vectors,  start  with  the  Af  —  1 
just  constructed  and  for  k  =  1  :  m  let  Qjk  =  0  and  let  j3jk  =  cvj  in  the  appropriate 
places.  We  illustrate  this  procedure  with  a  simple  example.  Consider  three  coarse 
strips  over  two  angles,  with  the  center  strip  of  each  view  refined  by  splitting  it  in  half. 
The  basis  is 

qi  =  (1  1  1  -  1  -  1  -  1  0  0  0  0)^ 

92  =  (10  1-1-1-1110  0)^ 

93  =  (111  -  1  0  -  1  0  0  -  1  -  1)^ 

Let  us  now  consider  the  matrix  5^^,  which  represent  only  the  fine  strip  func¬ 
tions.  This  region  is  depicted  for  four  views  in  the  figure  below.  Note  that  the  refined 
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Figure  55.  Refined  natural  pixels  for  four  angles. 


rays  do  not  cover  the  whole  image.  Only  the  refined  region  Qr  is  totally  covered  by 
fine  strip  functions  at  each  angle. 


Now,  for  a  vector  to  be  in  the  the  refined  natural  pixels  when 

overlayed  must  form  a  zero  image.  We  know  that  for  this  to  occur  in  fi/?,  must  be 
constant  by  angle  with  the  constants  summing  to  zero.  This  fact  provides  the  next 
result. 


Theorem  7.4:  is  of  full  rank. 

Proof:  Assume  there  exists  a  non-zero  vector  such  that  =  0. 

Then  v  must  be  constant  by  angle  with  oci  =  0.  If  all  the  fine  strips  are 
superimposed  over  the  image  at  once,  they  will  divide  it  into  a  collection  of 
polygons.  By  assumption,  only  the  refinement  region  Qr  is  completely  covered 
by  strips  at  each  angle.  Therefore,  there  must  exist  a  polygon  Pk  that  borders 
CIr,  and  in  particular  is  adjacent  to  some  polygon  P,  €  0,r,  but  is  not  formed 
from  any  of  the  strips  at  angle  (f>j.  Since  is  constant  by  angle  with  constants 
summing  to  zero,  then  for  these  two  polygons  we  must  have 

M-l  M-l 

ae  = 

e=i  e=i 
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Figure  56.  The  strip  functions  for  a  two-level  refinement. 

which  implies  that  aj  =  0.  Since  was  arbitrarily  chosen,  all  the  aj  must 
be  zero.  Therefore  =  0,  which  is  a  contradiction.  Hence  must  be  of 

full  rank.  | 

An  immediate  consequence  of  this  theorem,  when  combined  with  (7.3)  is 
Corollary  7.4:  If  then  G 

All  of  the  preceding  results  are  derived  from  a  spotlight  CT  problem  involv¬ 
ing  one  level  of  refinement.  .Most  will  generalize  to  multiple  refinement  levels.  For 
example,  assume  that  a  portion  of  the  original  refinement  region  is  itself  refined,  pro¬ 
ducing  two  levels  of  refinement.  The  strip  functions  for  two  angles  of  such  a  two  level 
refinement  are  shown  in  the  Figure  56. 

Letting  the  global  coarse  level  be  denoted  as  4/i,  the  first  refinement  level  as 
2h,  and  the  finest  level  as  h,  then  the  composite  system  matrix  will  have  the  form 
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Q4h4h 

Q4h2h 

Q4hh 

^2h4h 

Q2h2h 

Q2hh 

Qh4h 

Qh2h 

Qhh 

This  matrix  is  also  non-negative,  symmetric  and  positive  semi-definite.  Each  of  the 
nine  blocks  comprising  it  have  themselves  a  M  X  M  block  structure,  where  M  is  the 
number  of  angles  used  in  x-raying  the  image.  The  same  summability  properties  are 
present  as  well.  That  is,  the  level  h  rows  of  the  matrix  formed  from  refining  a  level  2h 
strip  function,  when  summed,  will  equal  that  row  corrsponding  to  the  level  2h  strip 
function  in  question.  Likewise,  level  2h  rows  will  sum  to  equal  the  level  Ah  rows  they 
are  refined  from.  This  allows  an  extension  of  Theorem  7.3,  which  we  believe  can  be 
proved  using  identical  arguments. 


Conjecture:  Let  B-  G  'JI^-^p+Q  the  composite  natural  pixel  discretized 
matrix  formed  at  M  angles,  and  assume  that  the  first  refinement  region  is 
covered  by  a  total  of  m\  level  4h  strip  functions  subdivided  into  level  2h  sti'ip 
functions,  and  assume  the  second  refinement  region  is  covered  by  a  total  of 
7712  level  2h  strip  functions  subdivided  mto  level  h  strip  fu7ictio7is,  as  outlmed 
above.  The7x  the  rank  of  B-  is  N  -{■  P  Q  —  (A/  +  mi  -|-  m2  —  1). 


This  recursive  refinement  can  be  expanded  to  as  many  levels  required,  in  the  same 
fzishion  as  the  second  level  was  added. 

From  our  analysis,  we  see  that  the  composite  matrix  B-  retains  many  of  the 
properties  of  the  unrefined  matrix.  Additionally,  we  have  determined  its  rank  and 
characterized  its  null  space.  In  the  next  section,  we  propose  a  multilevel  solution 
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technique  for  B-a-  =  /-  that  will  ultimately  be  equivalent  to  the  well-known  Fast 
Adaptive  Composite  (FAC)  method. 

C.  MULTILEVEL  APPROACH 

The  spotlight  CT  problem  is  a  composite  grid  problem,  in  which  an  operator 
equation  Lu  =  /  must  be  solved  on  some  composite  grid  17-  comprised  of  a  global 
coarse  grid  and  a  local  fine  grid  17^  (which  itself  could  be  a  composite  grid, 
allowing  for  recursive  refinement).  Fast  Adaptive  Composite  grid  methods  (FAC) 
were  developed  to  utilize  multilevel  technology  to  solve  such  composite  problems  in 
an  efficient  manner  [Ref.  34,  3]. 

FAC  methods  are  characterized  by  their  use  of  a  composite  grid,  which  is  the 
union  of  regular  grids  of  various  sizes.  The  problem  is  discretized  and  solved  on  the 
non-uniform  composite  grid,  but  all  of  the  actual  computations  occur  on  the  uniform 
subgrids.  This  provides  the  advantage  of  using  existing  uniform  grid  solvers,  while 
at  the  same  time  allowing  for  effective  resoltion  of  local  areas  of  interest.  For  this 
reason,  FAC  is  preferable  to  just  solving  the  system  of  composite  equations. 

There  are  three  main  features  that  allow  FAC  to  handle  grid  refinement  prob¬ 
lems  successfully.  First,  the  composite  grid  is  the  union  of  a  sequence  of  nested 
uniform  grids,  which  simplifies  the  data  structure  needed  to  represent  it.  Second, 
almost  all  computation  is  restricted  to  these  uniform  grids.  Finally,  the  use  of  multi¬ 
level  processing  to  correct  coarse  grid  approximations  with  fine  grid  residuals  through 
the  use  of  overlapping  grids  and  interpolation  at  the  grid  interfaces  allows  for  effective 
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intergrid  communication. 

As  with  all  multilevel  methods,  it  is  necessary  for  quantities  to  be  accurately 
represented  on  the  various  grids,  and  intergrid  transfer  operators  must  exist  to  transfer 
these  quantities  between  grids.  FAC  is  no  exception.  Letting  superscripts  denote  the 
grid  on  which  a  quantity  is  defined,  the  composite  grid  equation  becomes 

=  fK 

Likewise,  the  equations  and  denote  the  problem  restricted 

to  the  global  coarse  grid  and  the  local  fine  grid,  respectively.  Assume  that  intergrid 
transfer  operators  Ip  ;  fi-  — )•  :  fl-  — >  fi-,  and  fi- 

exist  to  transfer  quantities  between  grids.  The  details  involved  in  deriving  these 
representations  and  operators  may  be  very  cumbersome.  In-depth  treatments  can  be 
found  in  [Ref.  35,  34,  3].  Once  all  of  these  components  are  in  place,  FAC  is  given  by 
the  following  steps. 

•  Set  =  If{p  -  lAu^) 

•  Solve 

•  Correct  u-  u-+ 

•  Set  =  /^(/-  —  LAu-) 

•  Solve 

•  Correct  u-  <—  u-  + 

In  general,  FAC  first  solves  the  restriction  of  the  composite  residual  equation  to  the 
global  coarse  grid,  using  this  solution  to  correct  the  composite  grid  approximation. 
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The  restriction  of  the  composite  residual  equation  to  the  local  fine  grid  is  then  solved 
and  the  solution  is  used  to  correct  the  composite  grid  approximation.  Although 
formally  the  procedure  calls  for  exact  solvers  on  both  the  coarse  and  fine  grids,  in 
practice  iterative  methods  or  other  multilevel  solvers  are  used. 

We  now  show  that  FAC  applied  to  the  spotlight  CT  problem  is  formally  equiv¬ 
alent  to  a  block  Gauss-Seidel  formulation  of  the  same  problem.  Note  that  a  composite 
grid  element  u-  G  given  by 


u 


2h 


\ 


/ 


ran  be  decomposed  as 

^  (7.5) 

liitergrid  transfer  operators  that  satisfy  (7.5 1  arc  given  by 


\ 

\ 

Is 

i 

0 

= 

and  /.:  = 

K  0  j 

(  //•  / 

where  Is  and  Ip  are  identity  operators  of  tin*  appropriate  sizes.  It  is  significant  to 
note  that  FAC  generally  does  not  have  such  simple  intergrid  transfer  operators.  The 
simplicity  in  our  case  is  a  direct  result  of  the  discretization  by  natural  pixels,  and  the 
fact  that  we  require  refinement  to  be  an  exact  partition  of  coarse  rays. 

Now,  consider  the  FAC  scheme  applied  to  B-u-  =  /-.  Initially,  we  compute 
the  residual  of  the  composite  problem  restricted  to  the  global  coarse  grid  as 
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Q2h2h  Q2hh 
Qh2h  Qhh 


\  (  2/.  ^  ^ 


1,^  , 
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=  /2/1  _  Q2h2h^2h  _  Q2hh^h^ 


which  implies  that 


/f  =  (In  0). 


Next,  the  coarse  grid  residual  error  is  computed  as  Note  that 

this  is  a  formal  treatment,  as  52/12/1  jg  singular.  The  current  approximation  u-  is  then 


corrected  as 


^ 


which  after  applying  (7.5)  becomes 


(7.6) 


We  now  compute  the  residual  of  the  composite  problem  restricted  to  the  local  fine 
grid  as 


r*  = 


^  h 


j2h 


rh 

J  / 


Q2h2h  Q2hh 
Qh2h  Qhh 


y  u*  y 
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which  implies  that 


4  =  (0  Jf)- 

The  fine  grid  residual  error  is  computed  as  r^,  and  the  current 

approximation  u-  is  again  corrected  as 


which  after  applying  (7.5)  and  (7.6)  becomes 

(_  /i(„*  +  e‘)  +  4 („“  +  e“).  (7.7) 


Now,  we  show  that  the  block  Gauss-Seidel  formulation  of  B-u-  =  f-  results 
in  (7.7)  as  well.  The  spotlight  CT  problem  can  be  expressed  as 


Q2h2h 

Qh2h 


Q2hh 

Bhh 


“  / 


J2h 

fH 


\ 

/ 


which  when  solved  for  and  yields  the  one  swt*ep  block  Gauss-Seidel  scheme 

•  .Set  4r- 

•  Set  <r-  {B^^)-\f^  - 


This  scheme  can  be  rearranged  to  produce  the  desired  results  in  the  following  fashion. 
For  the  global  coarse  grid  we  have 


(7.8) 

(7.9) 

(7.10) 
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Proceeding  in  the  same  manner  with  the  local  fine  grid  yields 


(7.11) 

(7.12) 

(7.13) 

The  argument  is  completed  by  applying  the  relationship  (7.5)  to  the  expressions  in 
(7.10)  and  (7.13)  above,  yielding 

^  =  4(e2'‘  +  u2^)  +  /^(e'‘  +  u'‘). 

This  shows  that  FAC  is  formally  equivalent  to  block  Gauss-Seidel  on  the  spot¬ 
light  CT  problem. 

The  methods  are  formally  equivalent,  because  we  already  know  that  the  block 
tnalrix  jg  singular.  In  practice,  to  utilize  the  block  Gauss-Seidel  approach  we 

solve  each  block  system  in  turn  with  an  iterative  method,  or  a  multilevel  method  such 
a-s  PMLV. 

A  considerable  body  of  theoretical  results  exists  for  two-level  FAC  [Ref.  3]. 
This  theory  requires  that  the  variational  properties  be  satisfied,  that  quantities  be 
meeisured  using  the  energy  norm,  and  that  the  operator  B-  be  positive  definite.  Hence 
and  B^^  are  non-singular.  Under  these  assumptions,  there  exist  convergence 
factors  for  FAC,  given  in  terms  of  the  spectral  radii  of  combinations  of  and  B^^ 
with  the  intergrid  transfer  operators.  These  convergence  factors  exist  for  the  case 
when  an  exact  solver  is  used,  and  for  the  case  when  relaxtion  is  used  to  approximate 
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the  exact  solver.  The  latter  case  is  also  a  function  of  the  relaxation  scheme  being 
used. 

The  equivalence  of  FAC  and  the  block  Gauss-Seidel  scheme  would  allow  this 
theory  to  be  directly  applied,  except  for  the  fact  that  our  operator  B-  is  not  positive 
definite.  Even  for  such  problems,  FAC  theory  may  apply  in  certain  circumstances. 
McCormick  [Ref.  3]  states  that 

Theorem  7.5:  If  B-  is  positive  semidefinite  and 

NS{B^)  C  n  {I^NSiB^'^)),  (7.14) 

then  existing  FAC  convergence  theory  is  applicable. 

Unfortunately,  (7.14)  is  not  satisfied  for  the  spotlight  CT  problem,  as  B^^  is  of  full 
rank  and  we  know  that  B-  is  rank  deficient. 

Even  without  this  theory,  numerical  results  are  promising.  In  the  reconstruc¬ 
tions  that  follow,  a  geometry  of  32  detectors  over  20  angles  is  used  for  the  global 
coarse  grid.  The  refinement  region  is  located  in  the  center  of  the  image,  as  depicted 
in  Figure  57  below. 

Sixteen  coarse  strip  functions  at  each  angle  are  refined  by  splitting  them  in 
half.  The  composite  grid  image  is  then  reconstructed  by  using  the  spotlight  (FAC) 
method  developed  above. 

The  first  two  examples  are  reconstructions  of  the  two  brain  phantoms  used  in 
prior  numerical  experiments.  (See  Figures  45  and  46).  In  each  case,  the  global  coarse 
grid  representation  on  the  left  lacks  detail.  The  geometric  objects  within  the  brain 
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Rermement 

Region 


Figure  57.  Location  of  refinement  region. 

are  difficult  to  identify.  The  local  fine  grid  representation  on  the  right,  at  twice  the 
resolution,  provides  better  detail  about  the  center  region  of  each  brain. 


In  the  third  example,  the  real  value  of  spotlight  CT  is  demonstrated.  A  brain 
phantom  containing  a  small  tumor  is  depicted  in  Figure  60.  Note  that  the  tumor  is 


Figure  58.  Global  and  refined  reconstructions  -  example  1. 
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Figure  59.  Global  and  refined  reconstructions  -  example  2. 


sufficiently  small  that  it  is  not  expected  that  we  see  it  on  the  global  coarse  reconstruc¬ 
tion.  We  reconstruct  this  image  using  the  spotlight  technique,  and  the  reconstruction 
is  shown  in  Figure  61.  In  the  global  coarse  reconstruction  to  the  left,  the  tumor  is  not 
apparent.  However,  the  higher  resolution  spotlight  image  on  the  right  clearly  shows 
the  presence  of  the  tumor. 

The  natural  pixel  discretization  of  the  spotlight  CT  problem  allows  a  high 
resolution  reconstruction  of  a  portion  of  an  image  at  a  lower  cost  than  that  required 
if  the  entire  image  is  discretized  at  a  fine  resolution.  When  solved  in  block  Gauss- 
Seidel  form,  further  savings  may  be  realized  by  using  the  PMLV  routine  as  a  solver. 
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Figure  60.  Actual  image  -  example  3. 


Figure  61.  Global  and  refined  reconstructions 


example  3. 
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VIII. 


CONCLUSIONS 


A.  GOALS  OF  THE  RESEARCH 

The  primary  goal  of  our  research  was  to  introduce  multilevel  technology  into 
the  solution  of  this  problem.  Intermediate  objectives  included  the  development  of  the 
mathematics  of  natural  pixels  more  fully,  development  of  spotlight  CT,  and  building 
a  better  mousetrap  -  development  of  a  multilevel  algorithm  that  is  competitive  with 
state-of-the-art  image  reconstruction  methods.  We  achieved  all  of  these  goals  but  the 
last  one,  and  laid  some  foundations  for  future  research. 

B.  STANDARD  ART 

We  focus  our  efforts  on  the  algebraic  reconstruction  technique  (ART),  one  of 
several  solution  techniques  commonly  used  to  solve  this  problem.  For  a  preliminary 
foundation,  the  standard  ART  approach  is  examined,  in  which  the  problem  is  dis¬ 
cretized  using  square  pixels  and  the  resulting  linear  system  of  equations  is  solved  using 
the  method  of  Kaczmarz.  We  learn  why  convergence  of  Kaczmarz  stalls  after  several 
iterations  on  the  large  rectangular  system  produced  by  this  dicretization.  Using  the 
Singular  Value  Decompostion  (SVD)  of  the  system  matrix,  it  is  determined  that  the 
singular  value  spectrum  can  be  separated  into  three  bands  -  a  resolvable  region,  a  near 
null  space,  and  a  null  space.  Numerical  tests  show  that  solution  components  in  the 
resolvable  region  of  the  spectrum  can  be  recovered  during  the  iteration,  while  those 
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in  the  near  null  space  cannot.  It  is  also  learned  that  Kaczmarz  mixes  modes  during 
iteration,  i.e.  singular  value  modes  over  the  entire  spectrum  are  excited,  including 
those  modes  in  the  near  null  space,  which  cannot  be  eliminated  later. 

C.  NATURAL  PIXEL  DISCRETIZATION 

We  adopt  the  natural  pixel  discretization,  in  which  the  image  is  discretized 
into  strips  corresponding  to  the  x-ray  paths  through  it,  on  each  of  which  the  image 
is  assumed  constant.  The  resulting  system  matrix  is  symmetric  and  positive  semi- 
definite,  which  allows  for  a  wider  variety  of  relaxation  methods  that  can  be  used 
to  solve  the  linear  system.  It  is  also  in  general  smaller  than  the  matrix  produced 
through  the  square  pixel  discretization.  We  take  advantage  of  the  symmetry  and 
reduced  matrix  size  and  solve  the  problem  faster,  while  still  reconstructing  high- 
quality  images.  A  detailed  linear  algebraic  examination  of  the  matrix  is  conducted, 
producing  several  useful  results.  The  rank  of  the  matrix  is  completely  determined  by 
the  geometry  used  to  x-ray  the  image,  and  the  null  space  of  the  matrix  is  characterized 
by  vectors  with  easily  recognized  properties.  We  construct  a  basis  for  the  null  space 
of  a  general  matrix  show  that  images  corresponding  to  such  vectors  are  invisible  and 
thus  do  not  affect  the  quality  of  the  reconstruction.  While  the  idea  of  natural  pixels 
did  not  originate  in  this  work,  the  analysis  of  the  matrix  properties  goes  well  beyond 
anything  previously  done,  and  the  spectral  analysis  is  entirely  new. 
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D.  GAUSS-SEIDEL  ITERATION 


The  Gauss-Seidel  method  is  considered  for  solving  the  linear  system.  We 
show  that  when  measured  with  the  energy  norm  Gauss-Seidel  cannot  diverge  on  this 
problem,  and  with  the  correct  initial  starting  vector  must  converge  to  the  minimum 
norm  solution.  We  examine  this  issue  further,  and  learn  that  the  eigenvectors  of  the 
Gauss-Seidel  iteration  matrix  associated  with  eigenvalues  of  modulus  one  are  in  the 
null  space  of  the  system  matrix.  This  fact  explains  why  convergence  is  ensured  for  the 
semidefinite  system.  We  find  that  Gauss-Seidel  applied  to  this  problem  stalls  after  a 
few  iteration  sweeps.  Using  the  techniques  developed  earlier,  the  performance  of  the 
iteration  is  analyzed.  As  with  the  Kaczinarz  matrix.  th<*  natural  pixel  system  matrix 
is  found  to  have  a  spectrum  that  separates  into  a  resolvable  region,  near  null  space 
and  null  space.  Gauss-Seidel  exhibits  behavior  similar  to  Karzinarz  in  that  it  cannot 
resolve  components  in  the  near  null  spare,  and  it  mixes  inodes.  Such  analysis  ha^  not 
btxM)  done  before  in  this  setting,  ami  help^  understand  why  the  iteration  stalls. 

E.  MULTILEVEL  METHODS 

This  type  of  behavior  is  often  exhibititl  by  (Jauss-Seidel  when  applied  to  sys¬ 
tems  resulting  from  the  discretization  of  PDK's  as  well.  It  is  well-known  that  multi¬ 
level  methods  improve  the  performance  of  relaxation  methods  such  as  Gauss-Seidel 
in  the  PDE  setting  [Ref.  28].  We  seek  a  similar  increase  in  performance  in  the  image 
reconstruction  setting. 

We  formally  cast  the  natural  pixel  discretized  image  reconstruction  problem  in 
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a  multilevel  setting  using  the  Multilevel  Projection  Methodology  (PML).  The  natural 
pixel  discretization  is  shown  to  be  a  discretization  by  orthogonal  projections,  which 
subsequently  induces  interspace  transfer  operators,  a  coarse  grid  correction  scheme 
and  a  relaxation  scheme.  This  is  the  first  time  PML  is  formally  applied  to  a  prob¬ 
lem  not  arising  from  partial  differential  equations,  and  represents  an  expansion  in 
the  types  of  problems  that  can  be  approached  with  multilevel  techniques.  A  PML 
V-cycle  algorithm  is  developed,  and  its  convergence  properties  established.  Numeri¬ 
cal  results  show  that  PML  initially  converges  faster  than  Clauss-Seidel  alone  on  the 
problem,  and  then  it  stalls  as  well.  While  the  bctifr  nioux Irap.a  fully  competetive 
multilevel  algorithm,  did  not  emerge,  nevertheless  th<’  I’.ML  method  can  solve  the 
problem  cheaper  and  faster  than  either  (Iauss-.Sel<lel  ahnie  *»r  I  he  standard  Kaczmarz 
approach,  while  producing  reconstruction''  t»f  comparable  quality.  This  represents  an 
improvement,  although  not  a  revolution,  in  the  siat«‘-of-tht'-art  of  the  algebraic  image 
reconstruction  problem. 

F.  SPOTLIGHT  CT 

Finally,  we  consider  the  Spotlight  (*T  problem,  where  high  resolution  is  desired 
for  only  a  portion  of  the  image.  We  discretize  tin*  problem  using  natural  pixels  on 
multiple  levels  of  resolution,  which  has  not  previously  been  attempted.  The  resulting 
composite  grid  avoids  the  high  cost  of  discretizing  the  whole  problem  on  a  fine  level. 
We  hope  to  solve  the  resulting  composite  linear  system  efficiently.  An  analysis  of  the 
system  again  reveals  a  rich  collection  of  properties.  As  in  the  one  level  case,  the  rank 
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rank  of  the  system  matrix  is  a  function  of  the  geometry  used  to  generate  it.  The  null 
space  is  once  again  characterized  by  vectors  with  easily  recognizable  properties,  and 
these  vectors  represent  invisible  images.  We  cast  the  composite  linear  system  in  a 
block  form  that  can  be  solved  using  a  block  Gauss-Seidel  scheme,  and  this  approach 
is  formally  shown  to  be  equivalent  to  the  multilevel  Fast  Adaptive  Composite  (FAC) 
method.  This  formulation  allows  us  to  solve  the  problem  on  uniform  grids  using 
the  techniques  developed  earlier,  instead  of  having  to  solve  the  composite  system. 
Numerical  results  for  two  levels  yield  high  quality  reconstructions. 

G-  FUTURE  RESEARCH 

Directions  for  future  research  would  concentrate  in  the  area  of  Spotlight  CT. 
Because  this  area  is  so  new  and  unexplored,  we  believe  that  it  is  here  our  results 
have  the  most  promise  of  making  a  positive  contribution.  Currently,  we  construct  the 
system  matrix  in  a  piecemeal  faishion,  one  block  at  a  time.  This  process  could  be  au¬ 
tomated,  and  the  discretization  could  be  taken  to  three  or  more  levels.  Theoretically, 
using  this  approach  an  image  could  be  resolved  to  as  fine  a  level  as  desired. 

Concurrent  research  on  the  image  reconstruction  problem  paralleling  our  own 
work  concentrated  on  a  similar  discretization,  but  with  the  image  region  defined  as 
that  region  in  the  intersection  of  all  the  views,  not  restricted  to  the  square.  With 
this  approach,  it  is  claimed  that  the  resulting  system  matrix  can  be  represented  by 
a  small  fraction  of  its  elements  [Ref.  23].  This  discretization  could  be  incorporated 
into  the  spotlight  CT  problem  as  well. 
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Finally,  the  open  question  of  finding  a  competitive  algorithm  based  on  multi¬ 
level  technology  remains.  This  would  be  a  great  accomplishment,  and  warrants  con¬ 
tinued  investigation.  Some  possible  avenues  to  pursue  include  coarsening  the  problem 
by  angles,  or  a  combination  of  angles  and  detectors.  These  are  but  a  few  of  the  areas 
that  appear  ideal  for  additional  research. 
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APPENDIX.  DETAILS  ON  GENERATING 

THE  MATRIX  B 


The  entries  of  the  matrix  B  are  the  areas  of  intersection  of  the  strips  defined 
by  the  x-rays.  There  are  N  such  strips,  each  of  which  defines  a  row  in  the  matrix.  To 
determine  the  N  entries  for  the  row,  we  must  compute  the  areas  of  intersection  of 
the  strip  with  all  N  strips  in  turn.  We  assume  the  image  is  contained  in  the  unit 
square,  so  we  are  only  concerned  with  the  intersections  of  strips  that  lie  within  this 
square. 

The  calculations  proceed  as  follows.  A  coordinate  system  is  imposed  on  the 
image  square,  with  the  origin  at  the  center  of  the  square.  There  are  two  cases  to 
consider  -  the  strips  are  parallel  so  that  if  they  intersect  the  area  is  just  the  area  of 
the  entire  strip,  or  the  strips  are  not  parallel,  in  which  case  they  intersect  in  the  form 
of  a  parallelogram.  In  the  former  case,  the  area  of  the  strip  is  calculated  bcised  on  how 
it  intersects  the  square.  There  are  six  ways  this  can  occur,  i.e.  it  contains  a  corner  of 
the  square,  it  contains  two  corners,  it  intersects  opposite  sides  of  the  square,  etc. 

If  the  strips  are  not  parallel,  then  we  compute  the  coordinates  of  the  vertices 
of  the  parallelogram  that  is  formed  from  their  intersection.  Next,  the  number  of 
vertices  that  lie  within  the  image  square  is  determined.  If  the  number  is  zero  or  four, 
the  calculations  to  find  the  area  are  trivial.  If  not,  then  we  must  determine  which 
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vertices  are  in  the  square  and  whether  or  not  the  intersection  contains  the  corner  of 
the  square.  The  area  is  then  calculated  based  on  these  determinations.  The  number 
of  unique  ways  the  parallelogram  can  intersect  the  border  of  the  square  are  far  too 
numerous  to  list  here.  Symmetry  is  exploited  to  reduce  the  number  of  calculations 
by  half,  that  is,  the  area  of  strip  i  with  strip  j  is  the  same  as  the  area  of  intersection 
of  strip  j  with  strip  i,  so  the  calculation  need  only  be  performed  once, 

A  brief  outline  of  the  algorithm  follows: 


Input:  m  =  #  of  angles 

nl  =  #  of  detectors  per  angle 


Output:  B,  an  N  X  N  matrix,  where  N  =  m+nl 


initialize  B=0 
for  i=l:m 
for  j=l :nl 

index=(i-l)*nl+j 
for  k=l :m 
for  l=l:nl 

indexl=(k-l)*nl+l 

if  indexl  >=  index  {exploit  symmetry} 

if  i=k  {rays  are  parallel} 


compute  area 

else  {rays  intersect} 

compute  coordinates  of  vertices 
determine  which  vertices  are  in  square 
compute  area 
end{else} 

B(index , indexl)=B(indexl , index)=area 
end{if } 
end{f or} 
end{f or} 
end{f or} 
end{f or} 
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